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ABSTRACT 


Pseudospectral  Collocation  Methods  for  the  Direct 
Transcription  of  Optimal  Control  Problems 

by 

Jesse  A.  Pietz 

This  thesis  is  concerned  with  the  study  of  pseudospectral  discretizations  of  optimal 
control  problems  governed  by  ordinary  differential  equations  and  with  their  applica¬ 
tion  to  the  solution  of  the  International  Space  Station  (ISS)  momentum  dumping 
problem. 

Pseudospectral  methods  are  used  to  transcribe  a  given  optimal  control  problem 
into  a  nonlinear  programming  problem.  Adjoint  estimates  are  presented  and  analyzed 
that  provide  approximations  of  the  original  adjoint  variables  using  Lagrange  multi¬ 
pliers  corresponding  to  the  discretized  optimal  control  problem.  These  adjoint  esti¬ 
mations  are  derived  for  a  broad  class  of  pseudospectral  discretizations  and  generalize 
the  previously  known  adjoint  estimation  procedure  for  the  Legendre  pseudospectral 
discretization.  The  error  between  the  desired  solution  to  the  infinite  dimensional  opti¬ 
mal  control  problem  and  the  solution  computed  using  pseudospectral  collocation  and 
nonlinear  programming  is  estimated  for  linear-quadratic  optimal  control  problems. 
Numerical  results  are  given  for  both  linear-quadratic  and  nonlinear  optimal  control 
problems. 

The  Legendre  pseudospectral  method  is  applied  to  formulations  of  the  ISS  momen¬ 
tum  dumping  problem.  Computed  solutions  are  verified  through  simulations  using 
adaptive  higher  order  integration  of  the  system  dynamics. 
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Chapter  1 
Introduction 

Optimal  control  problems  (OCPs)  governed  by  ordinary  differential  equations  arise  in 
a  wide  range  of  applications.  One  particular  field  where  these  optimal  control  prob¬ 
lems  are  abundant  is  the  aerospace  industry.  Aerospace  engineers  have  been  solving 
optimal  control  problems  for  trajectory  optimization,  spacecraft  attitude  control,  jet 
thruster  control,  missile  guidance  and  many  other  applications  for  decades.  Methods 
for  obtaining  these  solutions  are  almost  as  copious  as  the  applications  themselves. 

A  traditional  approach  to  solving  OCPs  entails  forming  the  optimality  conditions 
directly,  using  the  calculus  of  variations  and  Pontryagin’s  maximum  principle  [42], 
and  then  solving  the  resulting  equations  to  obtain  the  solution  to  the  optimal  control 
problem.  This  is  known  as  the  indirect  approach  for  solving  OCPs.  The  references 
[9,  23,  41,  42,  45,  46]  present  just  a  small  sample  of  the  work  that  discusses  or  applies 
indirect  methods  for  the  solution  of  optimal  control  problems.  In  rare  cases  the 
solution  can  be  obtained  analytically  from  these  optimality  conditions,  but  in  general, 
approximation  methods  are  used  to  solve  the  problem  numerically.  The  optimality 
conditions  of  these  problems  generally  take  the  form  of  differential  algebraic  equations 
(DAEs)  or  boundary  value  problems  (BVPs).  The  approximate  solution  to  the  OCP 
can  be  obtained  by  using  a  BVP  solver.  Many  such  methods  exist.  Perhaps  the  most 
popular  methods  are  multiple  shooting  and  collocation.  The  reader  is  encouraged  to 
consult  [3]  for  more  information  on  these  and  other  numerical  methods  for  solvings 
BVPs. 

Alternatively,  one  can  discretize  the  governing  ODEs  and  the  integral  terms  in  the 
objective  functional  or  constraint  functions  and  thereby  replace  the  infinite  dimen- 
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sional  optimal  control  problem  by  a  large  nonlinear  programming  problem  (NLP). 
This  is  known  as  the  direct  or  direct  transcription  approach  for  solving  OCPs.  This 
approach  is  typically  easier  to  use,  especially  for  OCPs  with  state  equality  or  inequal¬ 
ity  constraints.  Direct  methods  have  been  used,  e.g.,  in  [6,  7,  16,  44], 

This  thesis  focuses  on  a  class  of  direct  transcription  methods  in  which  the  govern¬ 
ing  ODEs  are  discretized  using  pseudospectral  collocation  methods.  Such  methods 
have  attracted  attention  [15,  14,  17,  18,  20,  43]  because  of  their  alleged  superior 
approximation  properties  and,  in  the  case  of  Legendre  pseudospectral  method,  the 
availability  of  a  so-called  adjoint  map  or  estimate.  However,  most  of  the  existing 
work  in  this  area  is  numerical  with  incomplete,  informal  discussions  of  mathematical 
properties  of  pseudospectral  discretizations  for  optimal  control  problems. 

The  goals  of  this  thesis  are  to  improve  the  mathematical  understanding  of  pseu¬ 
dospectral  discretizations  for  optimal  control  problems  and  to  apply  these  methods 
to  the  solution  of  optimal  control  problems  with  significance  to  the  aerospace  com¬ 
munity.  In  particular,  we  provide  a  systematic  derivation  of  adjoint  estimates  for  all 
pseudospectral  discretizations  that  use  Gauss-Lobatto  points  and  we  present  rigorous 
results  on  the  error  between  the  solution  computed  using  pseudospectral  discretiza¬ 
tions  and  the  exact  solution  of  the  underlying  infinite  dimensional  OCP. 

Adjoint  estimates  provide  approximations  to  the  adjoint  variables  (also  known 
as  costate  variables)  corresponding  to  the  optimal  solution  of  the  OCP  in  terms  of 
the  Lagrange  multipliers  corresponding  to  the  NLP  derived  using  the  direct  tran¬ 
scription  method.  Such  approximations  are  important  for  error  analysis,  mesh  refine¬ 
ment  strategies,  and  real-time  optimization  using  the  method  of  neighboring  extrema. 
Among  the  few  results  on  adjoint  estimates  are  [17,  26,  27].  In  the  context  of  pseu¬ 
dospectral  discretizations,  only  [17]  have  provided  an  adjoint  estimation  procedure  for 
the  particular  case  of  Legendre  pseudospectral  discretizations.  This  thesis  provides 
a  systematic  derivation  of  adjoint  estimates  for  all  Gauss-Lobatto  pseudospectral 
discretizations,  which  as  a  special  case  includes  the  result  of  [17].  The  work  on  ad- 
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joint  estimation  provides  the  foundation  towards  a  rigorous  convergence  analysis  that 
provides  estimates  for  the  error  between  the  solution  of  the  infinite  dimensional  opti¬ 
mal  control  problem  and  associated  adjoint  as  well  as  the  solution  of  the  discretized 
optimal  control  problem  and  associated  Lagrange  multipliers.  Such  error  estimates 
are  not  available  in  the  existing  literature.  This  thesis  derives  error  estimates  for 
linear-quadratic  optimal  control  problems,  and  presents  numerical  results  for  both 
linear-quadratic  and  nonlinear  optimal  control  example  problems. 

In  the  second  part  of  this  thesis,  a  class  of  pseudospectral  direct  transcription 
methods  are  applied  to  a  series  of  optimal  control  problems  derived  from  the  International 
Space  Station  (ISS)  momentum  dumping  problem.  This  is  an  attitude  control  prob¬ 
lem  where  the  attitude  of  the  station  is  manipulated  by  a  controller  which  uses  control 
moment  gyroscopes  (CMGs).  The  issue  here  is  that  the  CMGs  have  a  maximum  mo¬ 
mentum  threshold  which  cannot  be  exceeded.  Doing  so  will  result  in  loss  of  control  of 
the  vehicle.  The  goal  is  to  find  a  control  trajectory  that  will  maneuver  the  attitude  of 
the  ISS  from  some  initial  state  to  some  final  state  with  minimal  total  momentum  on 
the  CMGs,  obeying  the  system  dynamics  and  never  exceeding  the  momentum  thresh¬ 
old  along  the  way.  What  makes  this  problem  difficult  is  the  severe  nonlinearity  of  the 
problem  and  the  possible  discrete  nature  of  the  controls.  Related  spacecraft  control 
problems  are  discussed  in  [1,  5,  8,  13,  36,  38,  43,  45].  This  thesis  includes  a  study  of 
the  numerical  solution  to  the  ISS  momentum  dumping  problem  which  demonstrates 
the  utility  of  pseudospectral  methods  for  the  direct  transcription  of  optimal  control 
problems. 

This  thesis  is  organized  as  follows.  Chapter  2  states  the  general  form  of  the 
optimal  control  problems  that  will  be  considered,  their  corresponding  optimality  con¬ 
ditions  and  provides  some  examples  problems  that  will  be  used  throughout  this  the¬ 
sis.  Chapter  3  states  the  optimality  conditions  of  the  discretized  OCP,  describes 
how  adjoint  estimates  are  obtained  and  explores  some  of  the  consequences  of  using 
pseudospectral  methods  in  the  direct  transcription  of  optimal  control  problems.  The 
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application  of  the  Legendre  pseudospectral  method  to  the  space  station  momentum 
dumping  problem  is  addressed  in  Chapter  4.  Finally,  Chapter  5  contains  remarks, 
conclusions  and  suggestions  for  future  work. 
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Chapter  2 


Optimal  Control  Problems  Governed  by  Ordinary 

Differential  Equations 


This  chapter  provides  a  description  of  the  optimal  control  problems  that  are  consid¬ 
ered  in  this  thesis.  In  addition,  first  order  necessary  optimality  conditions  are  stated 
for  the  infinite  dimensional  problem  and  reformulations  of  optimal  control  problems 
are  presented.  This  material  is  well  known  and  will  be  used  in  subsequent  chapters. 
Much  of  the  material  developed  in  Chapter  3  can  be  applied  to  OCPs  that  are  of  a 
more  general  form  than  (2.1).  However  (2.1)  covers  the  applications  considered  in  this 
thesis  and  is  sufficient  to  describe  the  theoretical  aspects  of  pseudospectral  methods 
as  used  to  transcribe  optimal  control  problems. 


2.1  Problem  Statement  and  Necessary  Optimality  Conditions 


In  this  thesis  we  consider  optimal  control  problems  of  the  following  class 

rtf 

min  m(y(tf))+  £(y(t),u(t))dt, 

Jt0 

s.t. 

jty{t)  =  f(y(t),u(t)), 

y(t0)  =  y0, 
b(y(tf))  =  o. 


(2.1a) 


(2.1b) 

(2.1c) 

(2-ld) 


Here  y  :  [to,tf]  h- >•  Mni/  are  the  state  variables  and  u  :  [t0,t/]  i— >•  are  the  control 
variables  to  be  determined.  The  functions  m  :  KC1'  i— »  M,  £  :  x  Mn“  i— »  M, 

/  :  x  i — *  IC*'  and  b  :  i— >•  Mnt  as  well  as  yo  G  are  given.  The  conditions 
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(2.1b)  and  (2.1c)  are  called  the  state  equations.  The  problem  (2.1)  is  said  to  be  in 
Bolza  form. 

We  seek  solutions  u  G  L^[t0,tf]  and  y  G  [to,  tf],  The  space  Loothn^/]  is 
the  space  of  all  Lebesgue-measurablc  functions  with  the  property  that  their  absolute 
value  is  essentially  bounded  on  [f0,£/]  and  L^[to,tf]  =  (Loo[t0,tf])nu ,  while 


wi,yooK  tf]  :=  {y  ■  [t0,tf] 


y  is  absolutely  continuous, 


These  spaces  are  equipped  with  norms 


MlliS?[to,t/]  =  ess  supIK^Ih 


and 

WyWw^tf]  =  max{||y|| iS[t0it/],  \\jty\\L^[t0,tf]}- 

The  necessary  optimality  conditions  for  (2.1)  can  be  obtained  using  a  general¬ 
ization  of  the  well-known  Lagrange  multiplier  theorem.  We  need  a  constraint  qual¬ 
ification  to  ensure  the  surjectivity  of  the  linearized  constraints  (2.1b)-(2.1d)  at  the 
solution  y* ,  u* .  For  (2.1)  this  can  be  guaranteed  by  assuming  that  the  Jacobian 
bv(y*(tf))  has  full  row  rank  and  that  the  linearized  state  equations 

jty{t)  =  fy(y*(t),u*(t))Ty(t)  +  fu(y*(t),u*(t))Tu(t), 
almost  everywhere  on  [t0,i/], 

y(to)  =  Vo,  (2-2) 


are  controllable.  Controllable  means  that  for  every  tjf  G  M.ny  there  exists  a  control 
u  G  L^s\to,tf]  such  that  the  solution  y  G  ,tf]  of  (2.2)  satisfies  y(tf)  =  y$. 

The  following  theorem  is  proven  in  [33,  Sec.  5.4], 

Theorem  2.1  Let  the  optimal  control  problem  (2.1)  be  given.  Let  the 
functions  m,  £,  f  and  b  be  continuously  partially  differentiable.  Let  w*  G 
L^[to,tf]  be  an  optimal  control  and  let  y*  G  ,tf]  be  the  resulting 

state.  Let  the  matrix  by(y*(tf))  have  full  row  rank  and  let  the  linearized 
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system  (2.2)  be  controllable.  Then  there  exists  a  function  p*  G  £/] 


and  vectors  (go)* 

G  Rny  and  (g/)*  G  Rnb  such  that  the  boundary  value 

problem 

d  .  . 

=  f{y*{t),u*{t)), 

y*(to) 

=  Vo, 

(2.3a) 

b(y*(tf )) 

=  o, 

d  .  . 

=  -fv(y*(t),u*(t))Tp*(t)  -  £v(y*(t),u*(t)), 

(2.3b) 

(adjoint  equation) 

P*(to ) 

=  -(go)*, 

P*(tf) 

=  mv{y*{tf))  +  by(y*(tf))T  (qf)*. 

(2.3c) 

(transversality  conditions) 

0 

=  fu(y*(t),u*(t))Tp*(t)  +  L(y*(t),u*(t)), 

(2.3d) 

(local  Pontryagin  maximum  principle) 

are  satisfied  almost  everywhere  on  [£0,t/]. 

The  conditions  (2.3b), (2.3c)  and  (2.3d)  are  obtained  by  computing  the  Frechet 
derivatives  of  the  Lagrangian  function 

£{y,u,p,q0,qf)  =  m(y(tf))  +  b(y(tf))Tqf  +  (y(t0)-y0)Tq0 

+  f  +  p(t)T[f(y(t),u(t))  -  jty(t)}  dt,  (2.4) 

with  respect  to  y  and  u  and  setting  these  Frechet  derivatives  to  zero. 


2.2  Bolza  and  Mayer  Forms  of  the  Optimal  Control  Problem 


For  most  of  this  thesis,  optimal  control  problems  of  the  following  Mayer  form  are 
considered, 


min  m(y(tf)),  (2.5a) 

s.t. 

=  f(y(t),u(t)),  (2.5b) 

y(to)  =  Vo,  (2.5c) 

b(y(tf))  =  0.  (2.5d) 


This  is  no  restriction,  since  every  problem  (2.1)  in  Bolza  form  can  be  converted  into 
an  equivalent  problem  in  Mayer  form.  This  will  be  discussed  shortly. 

The  optimality  conditions  for  the  Mayer  form  optimal  control  problem  can  be 
obtained  as  an  application  of  Theorem  2.1.  They  are  stated  here  for  later  reference. 
They  consist  of  the  boundary  value  problem 


=  f{y{t),u{t)), 

y(t  o) 

=  Vo , 

(2.6a) 

Ky(tf)) 

=  o, 

the  adjoint  equation 

=  -fyi.y(t),u(t))Tp(t),  (2.6b) 

with  transversality  conditions 

p(to)  =  -Qo, 

P(tf)  =  my(y{tf))  +by(y(tf))Tqf,  (2.6c) 

and  the  gradient  equation 


fu(y(t),u(t))Tp(t)  =  0. 


(2.6d) 
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One  can  transform  a  Bolza  problem  into  a  Mayer  problem  by  moving  the  integral 
term  in  (2.1)  into  the  differential  equation  by  defining  an  auxiliary  variable 

z(t)  =  It0  Ky(r),u(r))dT.  (2.7) 

If  (2.7)  is  differentiated  with  respect  to  t  the  following  initial  value  problem  (IVP)  for 
z  arises 


z{to)  =  0. 

(2.8) 

The  IVP  (2.7)  could  then  be  inserted  into  the  differential  equation  in 

(2.5b)  along 

with  y.  This  leads  to  the  following  optimal  control  problem  (2.9)  in 

which  is  equivalent  to  (2.1). 

Mayer  form, 

min  m(y(tf))  +  z(tf), 

(2.9a) 

s.t. 

jty{t)  =  f(y(t),u(t)), 

(2.9b) 

Jtz(t)  =  £(y(t),u(t)), 

(2.9c) 

y(to)  =  Vo, 

(2.9d) 

z(t0 )  =  0, 

(2.9e) 

Ky(tf))  =  o. 

(2.9f) 

Application  of  Theorem  2.1  to  (2.9)  leads  to  the  following  necessary  optimality 
conditions.  There  exist  adjoint  variables  p  associated  with  (2.9b)  and  r  associated 
with  (2.9c)  as  well  as  multipliers  q0,q f  associated  with  (2.9d),  (2.9f)  and  multipliers 
s0  associated  with  (2.9e)  such  that  the  constraints 

£y(t)  =  f(y(t),u(t)), 
y(t  o)  =  Vo, 

Ky(tf))  =  °> 


(2.10a) 


and 
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i  z(t)  =  *(y(t),u(t)), 

z(t0)  =  0, 

are  satisfied,  the  adjoint  equation 

jtP(t)  =  ~fy(y(t)^(t))Tp(t)  ^£y(y(t),u(t))r(t), 

the  auxiliary  adjoint  equation 

d  .  . 

dtr{t)  =  °- 

the  transversality  conditions 

p(to)  =  -go, 

p{tf)  =  my(y{tf))  +  by{y{tf))Tqf, 


(2.10b) 


(2.10c) 


(2.10c!) 


(2.10e) 


the  auxiliary  transversality  conditions 

r(«0)  =  -So, 

r(tf)  =  1,  (2.1Qf) 

are  satisfied,  and  the  gradient  equation 

fu(y(t),u(t))Tp(t)  +  £u(y(t),u(t))r(t)  =  0,  (2.10g) 

holds. 

Note  that  (2.10d)  and  (2.10f)  imply  r(t)  =  1.  The  necessary  optimality  conditions 
for  (2.1)  and  (2.9)  are  equivalent. 

The  conversion  of  the  problem  (2.1)  in  Bolza  form  into  a  problem  (2.9)  in  Mayer 
form  is  also  important  from  a  numerical  point  of  view.  For  an  efficient  numerical 
solution  it  is  important  that  the  discretization  of  the  state  equation  (2.1b)  is  consistent 
with  the  discretization  of  the  integral  in  (2.1a).  This  is  not  straightforward  for  many 
high  order  discretization  methods.  This  difficulty  is  avoided  when  (2.1)  is  transformed 
into  (2.9),  since  only  a  system  of  differential  equations  has  to  be  discretized.  The 
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discretization  of  the  ^-component  of  this  system  implicitly  defines  a  discretization  of 
the  integral  term  in  (2.1a)  that  is  consistent  with  the  discretization  of  (2.1b).  For 
additional  discussions  see,  e.g.,  [7]. 


2.3  Prototypical  Examples 

Throughout  this  thesis  two  example  problems  are  used  to  demonstrate  various  prop¬ 
erties  associated  with  solving  optimal  control  problems  using  a  pseudospectral  direct 
transcription  method.  These  problems  are  stated  here  so  that  they  may  be  referred 
to  elsewhere. 


Example  2.2  (Linear- Quadratic  Optimal  Control  Problem)  This  prob¬ 
lem  was  adapted  from  [27].  Consider  the  following  linear-quadratic  opti¬ 
mal  control  problem 


min  f*  y(t)2  +  \u(t)2  dt, 
s.t. 

(2.11) 

fty(t)  =  \y(t)pu(t ),  te  [0,1], 

j/(o)  =  i. 


Simple  evaluation  of  the  necessary  conditions  (2.3)  leads  to  the  following 
exact  solution  for  the  state 


oe3t  +  e3 

*<*>  =  em¥^ry  t6M' 

the  control 

2  (p3t  — 

=  e3«/2(2  +  e3)  ’ 

and  the  adjoint 


(2.12a) 


(2.12b) 


P*(t) 


2(e3t  -  e3) 
e3t/2(2  +  e3)  ’ 


te  [0,1]. 


(2.12c) 
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The  problem  (2.11)  can  be  equivalently  written  as  a  problem  in  Mayer 
from: 

min  1/2(1), 


s.t. 


ity i(*)  =  lyi(t)  +  u(t),  t  e  [0, 1], 

iMt)  =  yi(t)2  +  \u(t)2,  te  [0,1], 

2/i(0)  =  1, 

1/2(0)  =  0. 


(2.13) 


Evaluation  of  the  necessary  conditions  (2.6)  leads  to  the  following  exact 
solution  for  the  state 


(yi)*(t) 


2e3t  +  e3  .  1 

e3t/2(2  +  e3)’ 


(2.14a) 


the  auxiliary  state 


(2/2)*  (*) 


e~3t(2eet  +  (e6  -  2)e3t)  -  e6 
(2  +  e3)2 


t  e  [0,1], 


(2.14b) 


the  control 

=  ^[O,!],  (2.14c) 

the  adjoint 

(PMt)  =  -J/2{2  +  eiy  tel  0,1],  (2.14d) 

and  the  auxiliary  adjoint 


(P2).(*)  =  1,  te[0,l],  (2.14e) 


Example  2.3  (Orbit  Transfer  Optimal  Control  Problem)  This  example 
is  adapted  from  [8].  This  problem  is  frequently  used  in  the  context  of 
pseudospectral  direct  transcription  methods  for  optimal  control  problems, 
see  [15,  14,  IT].  Consider  the  following  orbit  transfer  optimal  control 
problem  of  finding  an  optimal  trajectory  and  thrust  steering  vector  to 


13 


transfer  a  spacecraft  from  an  initial  orbit  to  a  final  orbit  in  a  fixed  amount 
of  time.  This  problem  is  stated  as 

min  -^2(50)2  -  §y3(50)2  +  j/i(50)_1, 
s.t. 

jtVi (t)  =  V2 (t)i  t  G  [0,  50], 
iMt)  =  ^-^  +  0.01sin(«(t)), 

&Va(t)  =  -^«  +  0.01cos(n(t)), 

2/i  (0)  =  1.1, 

1/2(0)  =  0, 

2/3  (0)  =  1/vTI, 

where  y\  is  the  state  which  describes  radial  distance,  z/2  is  the  state  which 
describes  the  radial  component  of  velocity,  y 3  is  the  state  which  describes 
the  tangential  component  of  velocity,  and  u  is  the  controllable  thrust 
steering  angle.  It  should  be  noted  that  this  problem  has  no  analytical 
solution,  however  numerical  solutions  to  this  problem  can  be  found  in 
[15,  14,  17], 


t  e  [0,50], 
t  G  [0,50], 


(2.15) 
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Chapter  3 


Direct  Transcription  of  Optimal  Control  Problems 


In  this  chapter  we  describe  and  analyze  pseudospectral  collocation  discretization  of 
the  optimal  control  problem  (2.5).  Such  direct  transcription  methods,  especially 
Chebyshev  collocation  and  Legendre  collocation  methods  have  received  significant 
attention  recently  [14,  16,  17,  19,  29,  30].  One  reason  for  this  is  the  alleged  fast  con¬ 
vergence  of  the  solutions  of  discretized  optimal  control  problem  to  the  solution  of  the 
underlying  infinite  dimensional  control  problem.  Another  reason  for  the  large  interest 
in  direct  transcription  methods  based  on  Legendre  collocation  is  the  availability  of 
an  adjoint  estimate  that  relates  the  Lagrange  multipliers  of  the  discretized  problem 
to  the  adjoint  variables  p  evaluated  at  collocation  points. 

While  there  is  numerical  evidence  that  shows  fast  convergence  of  the  solutions 
of  discretized  optimal  control  problem  to  the  solution  of  the  underlying  infinite  di¬ 
mensional  control  problem,  the  theoretical  foundation  is  largely  missing.  The  papers 
[14,  17,  19]  cite  estimates  in  [10,  25]  for  errors  between  a  function  (its  derivatives) 
and  its  interpolant  (derivatives  of  its  interpolant).  Such  a  result  may  be  used  to 
establish  consistency  results  as  one  step  in  the  argument  that  pseudospectral  collo¬ 
cation  are  fast  converging  schemes  for  the  solution  of  the  dynamics  for  a  given  fixed 
control.  But  these  results  are  not  sufficient  to  establish  convergence  of  the  solution 
of  the  discretized  optimal  control  problem  to  the  solution  of  the  underlying  infinite 
dimensional  control  problem.  In  fact,  [27]  contains  simple  examples  using  Runge- 
Kutta  discretizations  of  optimal  control  problems,  which  show  that  the  solution  of 
the  discretized  optimal  control  problem  may  converge  to  the  solution  of  the  underly- 
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ing  infinite  dimensional  control  problem  with  a  slower  rate  than  one  might  expect,  or 
may  not  converge  at  all. 

The  goal  of  this  section  is  to  obtain  a  better  understanding  of  pseudospectral  col¬ 
location  discretizations  of  the  optimal  control  problem  (2.5),  in  particular  to  obtain 
a  better  theoretical  foundation  for  their  observed  fast  convergence.  An  important 
step  towards  this  goal  is  obtaining  estimates  for  the  error  between  the  true  solution 
to  the  OCP  and  the  computed  solution  to  the  OCP.  To  obtain  these  error  estimates, 
the  derivation  of  an  adjoint  estimate  that  relates  the  Lagrange  multipliers  of  the  dis¬ 
cretized  problem  to  the  adjoint  variables  p  evaluated  at  collocation  points  is  necessary. 
Such  adjoint  estimates  are  also  important  because  the  adjoint  variables  p  are  used, 
e.g.,  in  the  method  of  neighboring  extrema  for  real-time  optimal  control.  Despite  the 
theoretical  and  practical  importance  of  adjoint  estimates,  there  are  few  results  for 
high  order  discretizations.  We  will  derive,  in  a  systematic  way,  adjoint  estimates  for 
a  large  class  of  pseudospectral  collocation  discretizations  of  the  optimal  control  prob¬ 
lem  (2.5).  The  adjoint  estimation  procedure  of  [17]  is  a  special  case  of  our  treatment. 
We  also  derive  some  important  estimates  for  the  error  between  the  solution  of  the 
discretized  optimal  control  problem  and  the  solution  of  the  underlying  infinite  dimen¬ 
sional  control  problem  for  linear  quadratic  problems.  For  a  linear  quadratic  optimal 
control  problem  we  investigate  stability  results  numerically.  Finally,  we  comment  on 
the  discretization  of  optimal  control  problems  in  Bolza  form  using  pseudospectral 
collocation  discretizations. 

Since  the  reference  interval  [to,tf]  can  always  be  transformed  to  the  standard 
interval  [—1,1]  using  the  change  of  variables 

1 1— >■  —1  +  2 (t  —  —  to), 

it  is  assumed  that 

M/]  =  [-i,i] 

throughout  this  chapter,  except  for  Section  3.6. 


16 


3.1  Pseudospectral  Discretization  of  the  Optimal  Control  Problem 


For  the  numerical  solution  of  the  optimal  control  problem  we  discretize  (2.5)  using  a 
pseudospectral  collocation  method  with  collocation  points 

C0  =  — 1,  Ci,  .  .  .  ,  Cjv-1  G  (  —  1, 1),  Cat  =  1. 

The  state  y  is  approximated  by  the  polynomial 


where 


yN(t)  =  ^yMt), 

3=0 

(3,1) 

N  .  _ 

^■w=nc._c.’ 

i=i  3  1 

(3,2) 

j  =  0, . . . ,  N,  is  the  jth  Lagrange  polynomial.  Clearly, 

yN(ci)  =  yj ,  j  =  o,...,N. 


Furthermore, 


cj)  =  ^Zy^Mcj) 


ftyN(c  o) 


V  /  V  VN  ) 

where  D  =  D®Iny  G  Wny(N+1')xny(N+1'>  is  the  Kronecker  product  between  the  so-called 
differentiation  matrix  D  G  ]ghAr+1)x(Ar+i)  with  entries 

Djk=  ^fc(ci)j  j,k  =  l,...,N  +  l,  (3.4) 

and  the  ny  x  ny  identity  matrix  In  .  Recall  that  the  Kronecker  product  of  two  matrices 
A  G  Mmxn  and  B  G  Kix!  is  defined  as 


A®B  = 


AmB 


AlnB 


Ami  B  ■  ■  ■  A  vn,n.  B 
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To  discretize  the  state  equation  (2.5b),  we  substitute  y  by  yN  given  by  (3.1)  and 
require  that  the  ODE  (2.5b)  holds  at  the  collocation  points  Co, . . . ,  cat.  Furthermore, 
we  insert  (3.1)  into  the  objective  function  (2.5a),  the  initial  condition  (2.5c),  and  the 
final  time  condition  (2.5d).  This  leads  to  the  following  pseudospectral  collocation 
discretization  of  (2.5). 


mm 


m(yN), 


s.t. 


(  llu  ^ 


D 


\yN  J 

Vo  =  Vo, 

Kvn)  =  o. 


^  f(yo,u0)  \ 

V  f(VN,UN)  ) 


(3.5a) 


(3.5b) 


(3.5c) 

(3.5d) 


It  is  important  to  note  that,  following  [14,  17,  19]  and  others,  we  include  a  col¬ 
location  condition  at  Co  as  well  as  the  initial  condition  as  constraints  in  (3.5).  This 
is  different  from  other  direct  transcription  methods  based  on  collocation,  where  only 
the  collocation  conditions  at  Ci, . . . ,  Cjv  as  well  as  the  initial  condition  are  included  as 
constraints  (see  e.g.,  [4,  7,  44]).  This  also  seems  different  from  the  way  pseudospectral 
methods  are  used  to  discretize  boundary  value  problems,  where  one  also  eliminates 
collocation  conditions  at  c0,  Cjv,  depending  on  the  type  of  boundary  conditions  spec¬ 
ified  (see,  e.g.,  [21,  48,  49]). 

With  the  inclusion  of  the  collocation  condition  at  Co  it  is,  in  general,  not  possible 
to  solve  the  ny(N+2)  discretized  state  equations  (3.5b),  (3.5c)  for  the  ^(IV+l)  states 
y0l ,  jjn  given  controls  u0, . . .  ,un,  even  if  the  infinite  dimensional  state  equation 
(2.5b),  (2.5c)  has  a  unique  solution  y  for  given  control  u.  This  is  quite  different 
from  the  collocation  discretizations  [4,  7,  44],  where  the  discretized  state  equation 
consists  of  ny(N+l)  equations  and  where,  under  suitable  assumptions,  the  discretized 
state  equation  has  a  unique  solution  y0, . . . ,  y^,  given  controls  uq,  . . . ,  Un ■  Hence  the 
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discretization  (3.5b),  (3.5c)  of  the  state  equation  (2.5b),  (2.5c)  only  makes  sense  in 
the  context  of  optimal  control,  but  not  for  simulations. 

The  choice  of  including  the  collocation  condition  at  Cq  leads  to  nice  adjoint  esti¬ 
mation  properties,  which  will  be  discussed  in  Section  3.3.  On  the  other  hand,  it  is 
not  obvious  that  the  linearization  of  the  constraints  (3.5b)-(3.5d)  has  full  row  rank, 
even  if  the  constraint  qualification  in  Theorem  2.1  holds  for  the  infinite  dimensional 
problem.  It  is  possible  to  check  this  condition  a  priori  for  problems  with  linear  con¬ 
straints,  but  for  problems  with  nonlinear  constraints,  it  is  typical  to  assume  that  this 
constraint  qualification  holds  and  then  verify  it  after  computing  the  estimated  solu¬ 
tion,  see  [47].  In  the  next  example  we  examine  the  numerical  rank  of  the  linearized 
constraints  for  the  discretization  of  the  two  problems  stated  in  Examples  2.13  and 
2.15. 


Example  3.1  In  this  example,  the  numerical  rank  of  the  constraint 
Jacobian  in  (3.5),  evaluated  at  the  computed  solutions  and  using  Legendre 
pseudospectral  collocation,  are  investigated  for  the  problems  in  Examples 
2.13  and  2.15.  The  numerical  rank  is  investigated  by  inspecting  the  ratio 
between  smallest  singular  value  and  largest  singular  value  of  the  constraint 
Jacobian.  Figure  3.1  depicts  the  ratio  of  the  minimum  singular  value  of 
the  constraint  Jacobian  divided  by  the  maximum  singular  of  the  constraint 
Jacobian  value  for  different  N.  While  for  fixed  N  the  constraint  Jacobians 
have  full  rank,  Figure  3.1  indicates  that  one  should  expect  numerical  rank 
deficiency  for  large  N. 

We  conclude  this  section  by  stating  a  few  facts  about  the  two  pseudospectral  col¬ 
location  methods  that  have  been  used  for  the  direct  transcription  of  optimal  control 
problems  [14,  17,  19].  Details  about  the  computation  of  collocation  points,  cor¬ 
responding  differentiation  matrices  and  quadrature  weights  may  be  found,  e.g.,  in 
[10,  21,  34,  49], 
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Figure  3.1:  Ratio  Between  Minimum  Singular  Value  and  Maximum 
Singular  Value  of  the  Constraint  Jacobian  at  the  Solution  for 
Various  N.  -  Left:  Linear-Quadratic  OCP  in  Mayer  Form  - 
Right:  Orbit  Transfer  OCP 


Example  3.2  (Chebyshev-Gauss-Lobatto  Collocation)  The  Chebyshev- 
Gauss-Lobatto  collocation  points  are 


cj  ~~ 


(3.6) 


The  Cj’s  are  the  extrema  of  the  Chebyshev  polynomials  TN(t)  in  [—1, 1], 
The  Cj' s  are  also  the  roots  of  (1  —  f2)^T/v(f)  in  [—1, 1],  The  corresponding 
differentiation  matrix  is  given  by 


Dj  k  < 


V(-iy+fc 

( ck~cj ) 

j  ~f~  k, 

2N2+1 

6 

j  =  k  =  0, 

2N2+1 

6 

j  —  k  —  N, 

1  ck 

2  (1  -c2) 

otherwise, 

(3.7) 


where  £o  =  6v  =  2  and  £i  =  . . .  =  £n-i  =  L 

Note  that  the  points  c3  =  —  cos(jn/N)  are  defined  so  that  —  1  =  Co  < 
Ci  <  . . .  <  Cjv_ i  <  cp  —  1.  In  the  literature,  one  often  hnds  the  dehnition 


Cj  =  +  cos(j7r/p).  In  the  latter  case  the  signs  in  the  differentiation  matrix 
given  in  (3.7)  need  to  be  reversed. 

With  the  weighting  function 


20 


g(t)  =  1/Vl  -t2, 


and  the  quadrature  weights 


wj  =  /  g(t)^j(t)dt  = 


tt/N  1  <  7  <  TV  —  1 , 
n/(2N)  j  —  0,N, 


the  quadrature  formula 


g{t)h(t)dt  ~  E  Wjh(cj), 


is  exact  of  degree  2 N  —  1,  i.e., 


g(t)h(t)dt  =  Wjh(cj)  Vh  G  P2at-i([-1,  !])■ 


Example  3.3  (Legendre-Gauss-Lobatto  Collocation)  The  Legendre- Gauss- 
Lobatto  collocation  points  Cj,  j  =  0, . . . ,  N,  are  the  roots  of  (1  — 
where  L^{t)  is  the  Legendre  polynomial  of  degree  N .  The  differentiation 
matrix  for  the  Legendre-Gauss-Lobatto  collocation  is  given  by 


Djk 


LN(cj)  1 

-^iv(Cfc)  Cj—Ck 

-N(N+ 1) 

4 


N{N+l) 

4 


3  ~f~  k, 

j  =  k  =  0, 


j  —  k  —  N, 


(3.10) 


The  weighting  function  is 


0  otherwise. 


g(t)  =  i, 


and  the  quadrature  weights  are 


Wi  =  Jjj{t)dt=  n(n  +  1)[Ln{C])]2- 

The  quadrature  formula  (3.9)  is  exact  of  degree  2 N  —  1. 


(3.11) 


21 


3.2  Optimality  Conditions  for  the  Discretized  Optimal  Control 
Problem 


The  Lagrangian  corresponding  to  (3.5)  is  given  by 


where 


and 


L{y,  u,  A,  /x0,  hn)  =  m(yN) 


N 


N 


f  ( Uj  >  uj )  E  Dj^Vk 


+^o(vo  -  yo)  +  ^nKvn), 

y  =  (yo,---,yJr)T, 

U=  (Uq,...,Ux)T, 


(3.12) 


Differentiating  the  Lagrangian  (3.12)  with  respect  to  the  y^s  and  setting  the  deriva¬ 
tives  to  zero  gives  the  discrete  adjoint  equations 

fy(yo,  u0)T\0  —  J2k=o  Dk,o^k  =  — /d, 

fy^UjfXj-ZLoDkA  =  o,  j  =  1, N  —  1,  (3.13a) 

fy(VN,  Un)TXn  —  EfcLo  Uk'N^k  =  —by(yN)Ty,N  —  ‘my(yN)- 
Differentiating  the  Lagrangian  (3.12)  with  respect  to  the  u/ s  and  setting  the  deriva¬ 
tives  to  zero  gives  the  discrete  gradient  equations 

fu(Vj,  Uj)T\j  =  0,  j  =  0, . . . ,  N.  (3.13b) 


For  completeness  we  also  state  the  discrete  state  constraints 

f(yj,  ui )  -  EikLo  Dj,kyk  =  o,  j  =  o, . . . ,  n, 

yo-yo  =  o,  (3.13c) 

b(yN)  =  0, 
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which  are  obtained  by  differentiating  the  Lagrangian  (3.12)  with  respect  to  A j,  j  = 
0, . . . ,  N,  no  and  /xjv  and  setting  the  derivatives  to  zero. 

Note  that  the  discrete  adjoint  equations  (3.13a)  consist  only  of  n^iV+l)  equations 
for  the  ny(N  +  2)  +  rib  Lagrange  multipliers.  Therefore  the  discrete  adjoint  equations 
alone  do  not  specify  the  Lagrange  multipliers,  given  Uj,yj ,  j  =  0, ...  ,1V.  This  is  a 
consequence  of  the  fact  that  we  include  both  the  collocation  condition  at  Co  and  the 
initial  conditions  as  constraints  into  our  discretized  optimal  control  problem  (3.5)  and 
therefore  have  ny(N  +  2)  discrete  state  equations  (3.5b),  (3.5c)  for  ny(N  +  1)  discrete 
state  variables  yj,  j  —  0, . . . ,  N. 

3.3  Pseudospectral  Method  Adjoint  Estimation  Properties 

As  stated  earlier,  a  goal  of  this  thesis  is  to  derive  an  estimate  for  the  error  between 
the  solution  of  the  optimal  control  problem  (2.5)  and  the  solution  of  the  discretized 
problem  (3.5).  One  step  towards  this  goal  is  to  identify  a  suitable  adjoint  estimation 
procedure,  i.e.,  to  compare  the  Lagrange  multipliers  A  j  and  the  adjoints  p  evaluated 
at  the  collocation  points  cr  Such  adjoint  estimates  are  not  only  useful  for  the  error 
estimate  described  above.  For  many  applications  it  is  important  to  know  an  approx¬ 
imation  of  the  adjoint  variables  p  of  the  infinite  dimensional  control  problem  (2.5). 
Adjoint  information  is,  for  example,  important  to  adaptive  mesh  refinement  for  opti¬ 
mal  control  problems.  It  is  also  important  in  real-time  optimization  of  optimal  control 
problems  using  neighboring  extrema  [12,  35,  52],  Despite  its  importance,  there  are 
few  results  about  adjoint  estimation.  Adjoint  estimates  for  optimal  control  problems 
discretized  using  Runge-Kutta  methods  are  discussed  in  [27].  An  adjoint  estimation 
procedure  for  the  Legendre-Gauss-Lobatto  pseudospectral  discretization  of  optimal 
control  problems  is  presented  in  [17].  In  this  section,  we  will  derive  adjoint  estimates 
for  a  broad  class  of  pseudospectral  discretization  of  optimal  control  problems.  Our 
class  of  adjoint  estimates  includes  those  presented  in  [17]  as  a  special  case  and  is 
derived  in  a  systematic  way. 
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3.3.1  Integration  by  Parts  Approach 

Let  Co, . . . ,  Cat  be  the  collocation  points,  let  g  :  [—1, 1]  — >  M  be  a  positive  weighting 
function  and  let  w0, ...  ,wn  be  positive  weights  such  that  the  quadrature  formula 

N 


L 


g(t)h(t)dt  ~  E  Wjh{cj), 

3=0 


is  exact  of  degree  2 N  —  1,  i.e., 

ri  N 


i 


g(t)h(t)dt  =  X  Wjh(cj)  X/h  e  V2N- 1([— 1,  !])• 

3=0 


Instead  of  (3.12)  consider  the  weighted  Lagrangian 


(3.14) 


(3.15) 


Lw(y,  u,A,  ij,o,/in)  =  m(yN)  +  (y„  —  ya)TIMi  +  b(yN)TyN 

N  N 

+  X  /(%’>%)  -  X  DikVk 

3=0 


k= 0 


(3.16) 


Clearly,  if  A j  =  WjXj,  then  L(y,  u,  A,  no,  Hn)  =  Lw(y,  u,  A,  /i0,  /qv).  We  define  the 
polynomials 

N 

yN(t)  = 


k= 0 
N 


uN(t )  =  YukMt), 


and 


Note  that  since 


k= o 

N 


A N(t)  =  X 


k= 0 


N 


dyN(t )  =  Yy^Mt)  e  Pjv-i([-i,i]), 


dta  v  /  ^™dt 

k= 0 

and  A N(t)  G  Vn([~  1, 1]),  the  equation  (3.15)  implies  that 

N  N  N 


d 


X  ^  X  D^k  =  Xwi(AJV(ci))r^JV(ci) 


j=0  fc=0 


3=0 

i: 


dt' 

g(t)(XN(t))T^tyN(t)dt. 
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Furthermore,  (3.14)  yields  the  approximation 

N  N 

J2w3xJf(yj,ui)  =  ^2wj{*N(cj))Tf(yN(ci)>uN(cj)) 

3=0  3=0 

~  J  g(t){\N(t))Tf(yN(t),uN(t))dt. 

Hence,  one  may  interpret 


Lw(y,  u,  -A,  ho,  hJv) 


m[yN  (1))  +  (iT(-l)  -sb)J^  +  6(ir(i)r^ 

+  J  9(t)(\N{t))T  jtyN{i)  ~  f(yN(t),uN(t))  dt.  (3.17) 

This  suggests  the  following  relation  between  the  Lagrange  multipliers  X3  and  the 
weighted  Lagrange  multipliers  X3  of  the  direct  transcription  and  the  adjoint  variables 
P 


— A,  =  A,  « 


Wi 


P(cj) 
9{cj)  ' 


(3.18) 


The  relation  (3.18)  means  that  A 3/w3  =  X3  should  be  identified  with  p(cj) / g(cj).  In 
general  A  j/wj  =  X3  is  not  identical  to  p(cj)  /  g(cj).  However,  in  Section  3.4  we  will 
give  conditions  that  guarantee 


J_T  .  P(ci) 

w3  3  g(cj )  ’ 

as  N  — >  oo. 

The  necessary  optimality  conditions  associated  with  (3.16),  obtained  by  differen¬ 
tiating  the  weighted  Lagrangian  with  respect  to  each  variable  and  setting  the  result 
to  zero  are  given  by  the  weighted- discrete  adjoint  equations 

fy(yo,  U0)TW0X0  —  Ylk=0  DkoWkXk  =  — ho, 

fy(Vj,  UjfwjXj  -  Ej=o  DkjWkXk  =  o,  j  =  1, . . . ,  N  -  1,  (3.19a) 

fy{yN-,uN)TwNXN  —  J2k= o  ^kN^kXk  =  —  by(y]y)T  Pn  —  rny(yN), 
by  the  weighted- discrete  gradient  equations 


fu(Vj ,  UjfwjXj  =  0,  j  =  0, . . . ,  N, 


(3.19b) 
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and  by  the  weighted- discrete  boundary  value  problem 


wj  ( )  -  ££Lo  A’fcJ/fc)  =  0,  j  = 


Vo -Vo  =  0, 


(3.19c) 


b(yN)  =  0. 

Note  that  these  optimality  conditions  (3.19)  can  also  be  obtained  from  (3.13)  by 
substituting  Xj  =  WjXj  in  (3.13). 

The  following  notation  will  be  useful.  For  a  given  function  /  :  [—1, 1]  ->  1  we 
define  its  interpolating  polynomial 


pN(f)(t)  =J2f(cMt)- 


(3.20) 


Note  that 


PN(f)(ci)  =  f{ci),  i  —  0, . . .  ,N, 


-PN(f)(ci)  =  Y,Dijf(cj)>  i  —  0, . . . ,  N. 


(3.21) 


The  following  lemma  provides  a  discrete  integration  by  parts  formula. 


Lemma  3.1  Let  z0, . . . ,  be  arbitrary  and  dehne  zN(t )  =  z^t). 

For  any  continuously  differentiable  function  p  :  [—1,1]  — >  M  the  equation 

N  N 

i= 0  ’  j= 0 

N  N 

=  p(cn)zn  -  p{c0)z0  -  ^  ~J±)Zi  X]  D^p(cj) 

i= 0  j= 0 

,  t'  (t>  ,Pv*  pW\ 


+L^{Ph^~W) 


-z  z  d 


\pN{-){t)  JtpN(P)m 


r  i  ~N  t 


/  (t)  d 

+ J  i9(t)^-{pN(p)(t)-p(t))dt 


(3.22) 
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holds.  Furthermore, 


P(cN)SkN  —  P(c0)Sk0 

N 


— rH')2DkiP(.Ci)  +  Ch(N,p,g),  (3.23) 

»<c‘)  ts 

k  —  0, . . . ,  N,  where  Sjk  is  the  Kronecker  delta  and 

Ck(N,p,g)  =  I  g(t)  ^PN(p(t)  -  ^k(t)dt 

+j_m  ( PMjKt )  - ||)  f/Nm)dt 

+  f  (Pn (?)(*)  ~P(t))dt-  (3-24) 


Proof  The  definition  of  Pn  and  equation  (3.15)  imply 

N  ,  N  N 
\  P{Ci)  \  n 

Z^Wia(c  ) 
i=0  9^l)  7=0 

=  J  ^9  {t)PN(^){t)jtzN  (t)dt 

=  Lp(t)  itzN{t)dt  +  L3{t)  -  y§) 

f1  d 

=  p(cN)zN-p(c0)z0-  —p(t)zN(t)dt 

+  /.9W  (^(f)W-f|) 

Similarly, 

TV  TV 

t=o  7=0 

fL  zN  d 

=  /  g(t)PN(  —  )(t)-PN(p)(t)dt 

.7-1  .9  dt 

r1  d  f1  /  ?Ar('/'i\  d 

=  Id  W  Jt^dt  +  L^  -  iw)  dtPNij,)(t)it 

f  ^  (f\  (j 

+  J  i9(t)-^Tt(PN(p)(t)-p(t))dt. 


These  identities  imply  (3.22). 
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Equation  (3.23)  follows  from  the  choice  Zj  =  Sj^. 


□ 


Remark  3.4  If  g  —  1,  then 


,zN,  zN 


Pn(-)-  —  =  Pn(zn)-z*=  0 
9  9 

for  all  z0, . . . ,  Zn- 

The  properties  of  the  Lagrange  polynomials  imply 


N 


p  /V’fcw.x  Mci )  ,  ,.s  Mt) 

PK(7)(i)=S^y*(t)  = 


g(ck)' 


Remark  3.5  The  integrals  in  (3.22)  and  (3.24)  are  well  dehned  for 
g(t)  =  1  and  g(t)  =  1  / \/ 1  —  t2.  This  is  obvious  for  g(t)  =  1.  In  the  other 
case  it  can  be  seen  from  the  following  argument.  Let  g(t)  =  l/\/l  —  t2 
and  let  h  :  [—1, 1]  — >  M  be  a  continuously  differentiable  function  with 
/r(±  1)  =  0.  By  the  LHospital  rule, 


mn  h{t)g{t) 


lim  44L 

1 /g(t) 

Hm  W  - 

-g'(t)/g2(t) 


lim  /r'(t) 

t— »±i 


0. 


Hence  the  integrands  in  (3.22)  and  (3.24)  are  bounded  on  [—1, 1]  and 
continuous  on  (—1, 1). 


The  next  lemma  shows  that  the  adjoint  variable  p  divided  by  the  weighting  func¬ 
tion  g  satisfies  the  weighted-discrete  adjoint  equations  (3.19a)  with  an  error  that  is 
dependent  on  the  true  adjoint  p  and  on  N. 


Lemma  3.2  If  p  satisfies  the  adjoint  equation  (2.6b),  and  q0  and  qj 
satisfy  the  transversality  conditions  (2.6c),  then 


T...  P(C o)  V-  n  ...  P(Ck) 


fy(y(co),u(c0))  W0 


g{c0) 


Y  Dko™kj7 Y  =  -qo  +  Tq(N,p,  g), 


k= 0 
N 


fy(y(cj),  u{Cj))Tw~^  -  Y  Dkjwk^\  =  rj(N,p,  g ), 


k= 0 


fy(y(cN),u(cN))TwNj^. y  -  Y  DkNWkj0)  =  ~by(y(CN^T qf  -  my(yM) 


g(ck ) 


j  =  1, . . . ,  IV  —  1, 


k= 0 


+vn(N,p,  g), 


(3.25) 


where 


r“(N,p,g)  = 


Wj 


g(ci 


d 


d 


i, :pN(p)(cj )  -  -rAcj) 


dt 


dt 


+  ejA  1P1  g)i  (3-26) 


with  ej(N,p,g)  defined  in  Lemma  3.1. 


Proof  Use  equation  (3.23)  and  the  fact  that  p  satisfies  the  adjoint  equations  (2.6b) 
and  transversality  conditions  (2.6c)  to  deduce 


f  (  (  \  (  \\T  P(Cj )  n  P(ck) 

fy  {y\Cj)iU\Cj))  Wj  g(c  \  2-*/  DkjWk  n(r  \ 


k= 0 


gA) 


Wj 


g(cj 


N 


fy(y(cj)  i  u(cj))T  p(cj)  +  YDjkP(ck) 


k= 0 


-  p(cN)SjN  +  p(c0)Sj0  +  €j(N,p,g ) 


w 


g(c: 


J  '  TtPN^Cj^  -  ^~Aci) I  +  iutt  I  ~rAci)  +  fy(cP  y(cj)iu(cj))Tp(cj) 


d 

dt1 


Wj 


g(ci 


d 

dt1 


-p(cN)SjN  +  p(c0)Sj0  +  €j(N,p,  g) 


Wj 


d 


d 


gA 

+Qo  Sj0 


dtp"M(-Ci)  ~  dtp{Ci) 


-  my(y(l))T  +  by(y(l))Tqf  djN  +  ej(N,p,g ), 


for  j  =  0, . . . ,  N. 


□ 
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The  residual  terms  (3.26)  can  be  estimated  using  results  from  [10,  11].  To  state 
these  results  we  need  the  norms 

"/"»  =  (3.27) 

and 

II/IIL  =  tll^/lll  (3-28) 

k=0 

Theorem  3.6  i.  Let  c0,  ...,Cjv  and  w0,...,Wn  be  the  Chebyshev- 
Gauss-Lobatto  collocation  points  and  corresponding  quadrature  weights 
defined  in  Example  3.2  and  g{t )  =  1/a/1  —  t2.  If  /  :  [—1, 1]  — >  M  is  s- times 
continuously  differentiable,  then  there  exists  a  constant  C  independent  of 
/  and  N  such  that 


\PN(f)-f\Wg<CN~S\\f\U 

(3.29) 

PN(f)-f\\i,9<CN2-s\\f\\s,gi 

(3.30) 

and 

(N  2\ 

J  <CN*- 1|/||», (3.31) 

ii.  Let  Co, . . . ,  cn  and  wq,  . . . ,  wn  be  the  Legendre-Gauss-Lobatto  colloca¬ 
tion  points  and  corresponding  quadrature  weights  defined  in  Example  3.3 
and  g(t)  =  1.  If  /  :  [—1, 1]  — >  M  is  s- times  continuously  differentiable, 
then  there  exists  a  constant  C  independent  of  /  and  N  such  that 


^(/)-/llo,i<CiV1'2-||/||,,1, 

(3.32) 

^(Z)  - /lli.i  <  cjv5/2-*||/|U,, 

(3.33) 

and 


PN{f)(Cj) 


1/2 


<CN5/2~s\\f\\Stl 


(3.34) 
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Proof  Estimates  (3.29)-(3.31)  can  be  found  in  [10,  p.298].  Estimates  (3.32)-(3.34) 
can  be  found  in  [10,  pp. 293/294].  □ 


Corollary  3.1  i.  Let  Co,  ...,cjv  be  the  Chebyshev-Gauss-Lobatto  col¬ 
location  points  defined  in  Example  3.2  and  g(t )  =  \j\J\  —  t2.  If  /  : 
[—1, 1]  — >  M  is  .s'-times  continuously  differentiable,  s  >  2,  and  cr  >  0  then 

there  exists  a  constant  C  independent  of  /  and  N  such  that 
i v  \  1/2 

E(r“(Jv./.9))2 


o=o 
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+CN~C 
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(3.35) 
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ii.  Let  Co, . . . ,  cat  be  the  Legendre-Gauss-Lobatto  collocation  points  de¬ 
fined  in  Example  3.3  and  g(t )  =  1.  If  /  :  [—1, 1]  — »  M  is  s-times  continu¬ 
ously  differentiable,  s  >  5/2,  then  there  exists  a  constant  C  independent 

of  /  and  N  such  that 

tv  \  V2 

EKOV,/,1))2)  2  CA^-U/IL,.  (3.36) 
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Proof  In  this  proof  C  >  0  denotes  a  generic  constraint  independent  of  N  and  /. 
The  Cauchy-Schwarz  inequality  yields 
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There  exists  aC>0  independent  of  N  such  that  the  inverse  estimate 

ll^V’jllo.s  5;  CN  \\ipj\\ o,9 

holds  (see  equations  (9.4.4)  and  (9.5.4)  in  [10]).  Furthermore,  there  exists  a  constant 
C  >  0  such  that 


N 


1/2 


j  1 1 0,3  E  C  |  'y  '  WjlJjj  ( Cj )  —  Cy/vJj 


J=0 


(see  equation  (9.3.2)  in  [10]).  Finally, 

Vi 


9 


—  1 1 1  / 9 1 1  oo  |  iftj  1 1 0,g  —  |  V’j  I  0,3  E  Cy/Wj . 
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Using  Wj  G  (0, 1)  and  ]>U=0  wt  =  1,  the  previous  inequalities  imply  the  existence  of  a 
constant  C  independent  of  N  such  that 


N 


£ML<cv\  E 

3=0  j= 0 

Consequently,  there  exists  C  >  0  with 
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If  9  =  1,  then  Pjv {'ipj/g)  ~  ^ j/g  =  Pn(i/>j)  ~  i>j  =  0,  j  =  0, . . . ,  N. 

The  inequality  \\PN(f)\\i)g  <  \\f\\i,g  +  ||-Pjv(/)  -  f\\i,g  and  (3.30),  (3.33)  imply 

||Piv(/)||i,3<C'||/||s,3  ViV. 

Using  Wj/g2(cj )  <  1,  we  find  that 
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(3.38) 
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The  desired  estimates  now  follow  from  (3.37),  (3.38)  and  Theorem  3.6.  □ 

The  following  consistency  result  is  an  immediate  consequence  of  Lemma  3.2. 


Lemma  3.3  (Consistency)  Let  p  satisfy  the  adjoint  equation  (2.6b),  let 


qo  and  qj  satisfy  the  transversality  conditions  (2.6c),  and  let  A j,  j  = 
0,...,iV,  Ho,I^n  satisfy  (3.19a).  If  fy{y{cj),  u(cj))  =  fv{yj,Uj),  j  = 


0, . . . ,  iV,  by(y(cN ))  =  by(yN)  and  mv(y(cN))  =  mv(yN),  then 

fy(y(co),u(c0))TW0  -  A0^  -  Dk0Wk  ~ 

=  ~(go  ~Mo)  +rg(N,p,g), 

fy(y(cj)iu(cj))TWj  -  Ai)  -  DCWk  |  _  Afc) 

=  rj(N,p,  g)  j  =  1,...,N  -1, 


fy(y(cN),u(cN))TWN  ~  ^n)  -  ^  DkNWk 


k= 0 


=  -by(y(cN))T (qf  -  Hn)  -  my(y(cN))  +  r%(N,p,  g) 


where  r“(lV,p,  g),  j  —  0, . . . ,  N,  is  defined  in  (3.26). 


Proof  This  result  follows  immediately  by  subtracting  the  weighted  discrete  adjoint 
equations  (3.19a)  from  (3.25).  □ 


Note  that  since  the  discretized  optimal  control  problem  (3.5)  has  ny(N  +  2)  + 
constraints,  but  only  n.y(N  +  1)  state  variables,  there  are  only  ny(N  +  1)  discrete 
adjoint  equations  for  the  ny(N  +  2)  +  n b  Lagrange  multipliers  Ao,  •  ■ . ,  \n,  Po ,P>n- 
Hence,  the  Lagrange  multipliers  cannot  be  computed  from  (3.19a)  alone.  Therefore, 
it  not  possible  to  use  Lemma  3.3  and  a  stability  result  to  obtain  an  estimate  for  the 
error  between  p(cj)/g(cj )  and  A j/wj.  Such  an  estimate  will  be  obtained  in  Section 
3.4,  where  the  entire  optimality  system  is  considered. 

The  following  example  illustrates  the  adjoint  estimate  (3.18)  applied  to  the  orbit 


transfer  problem,  Example  2.3. 
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Example  3.7  Consider  the  Example  2.3.  We  apply  a  Legendre  and  a 
Chebyshev  psendospectral  discretization,  with  N  =  100,  to  this  problem. 
Figure  3.2  shows  the  Lagrange  multipliers  A j  as  well  as  estimated  adjoint 
variables  p(cj)  ~  g(cj)\j/wj  for  each  discretization.  Since  the  weighting 
function  g(t)  =  1/x/l  —  t2  for  the  Chebyshev  pseudospectral  methods  is 
singular  at  ±1,  the  estimated  adjoint  g(cj)\j/wj  becomes  less  accurate  as 
t  ->  ±1. 


(t)*g(t)  Lagrange  Multipliers 
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Figure  3.2:  Lagrange  Multipliers  XN  and  Estimated  Adjoints  gXN .  Top 
Left:  Legendre  Pseudospectral  Lagrange  Multipliers  XN . 
Top  Right:  Chebyshev  Pseudospectral  Lagrange  Multipliers 
XN .  Mid  Left:  Legendre  Pseudospectral  Adjoint  Estimates. 
Mid  Right:  Chebyshev  Pseudospectral  Adjoint  Estimates. 
Bottom  Middle:  Error  Between  Legendre  Adjoint  Estimates 
and  Chebyshev  Adjoint  Estimates 
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3.3.2  Weighting  Matrix  Minimization  Approach 


In  the  previous  section,  we  have  obtained  the  consistency  result  in  Lemma  3.3  by 
rewriting  the  adjoint  equations  (2.6b)  evaluated  at  Cj,  j  =  0, . . . ,  N,  in  the  form  of 
the  weighted  discrete  adjoint  equations  (3.19a)  using  the  discrete  integration  by  parts 
formula  (3.22). 

Alternatively,  one  may  consider  an  approach  which  is  motivated  by  the  adjoint 
estimation  procedure  in  [17].  In  this  approach,  we  identify 


j  =  0, . . . ,  A, 


w 


(3.39) 


where  Wj,  j  =  0 , ,N  are  suitably  chosen  weights  to  be  determined  below. 
Let  Wj  7^  0,  j  —  0, . . . ,  N,  and  consider  the  identity 


N 


fy  {Vj )  uj )  At  ^  ^  Dkj^k 


k= 0 


N 


fyil/j.  //,)'  A,  -  4 


k= 0 
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A/c 

wk 


N 


fy{Vjiuj)  \  Y'XWkPkj  +  wjDjk )  ~  +  ^  ]  WjDjk  ~ 

A:=0  Wk  k= 0  Wk 


=  VJi 


N 


fy(yj,uj)T^-  +  YD^Xk 


k= 0 


Wk 


N 


Yi^Dkj  +  WjDjk (3.40) 


k= 0 


Wk 


If 


'MkDfzj  +  WjDjk  ^ 


(3.41) 


-1  j  =  k  =  0, 

1  j  =  k  =  N, 

0  otherwise, 

then  we  will  show  below  that  A j/wj  and  p(cj )  are  related.  However,  the  identities 
(3.41)  cannot  always  be  satisfied.  Therefore,  let  E  €  ]gAAr+1)x(lv+1)  be  the  matrix  with 
entries 

'  -1  j  =  k  =  0, 

1  j  =  k  =  N, 

0  otherwise 


Ejk  =  \ 


(3.42) 
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and  choose  the  w/ s  such  that 

N 

(' wkDkj  +  w,Djk  -  Ejky  =  || WD  +  DtW  -  E\\% 

j,k=0 

is  minimized,  where  W  =  diag(w0, . . .  ,wn)  and 
problem 

min  ||  WD  +  DtW 

Wo,...,WN 

is  a  linear  least  squares  problem. 

Remark  3.8  For  any  choice  of  collocation  points  for  which  the  cor¬ 
responding  differentiation  matrix  D  has  a  nonzero  diagonal  entry  DJ3 , 
j  G  {1, . . . ,  N  —  1},  there  is  no  Wj  y  0  such  that 

wjDjj  =  ~wjDjj.  (3.44) 

Consequently,  in  this  case  there  are  no  Wj  y  0,  j  =  1, . . . ,  N  —  1,  for  which 
WD  +  DtW  -E  =  0. 

The  differentiation  matrix  (3.7)  for  the  Chebyshev  pseudospectral  method, 
satisfies  DJ3  y  0,  j  —  0, . . . ,  N. 

The  differentiation  matrix  for  the  Legendre  pseudospectral  method,  which  uses 
Legendre-Gauss-Lobatto  points,  satisfies  Du  =  ...  =  -D/v-qjv-i  =  0  (see  Example 
3.3).  In  this  case  the  solution  of  the  linear  least  squares  problem  (3.43)  is  known  and 
satisfies  (3.41). 

Lemma  3.4  Let  c3,  j  =  0, . . . ,  N,  be  the  Legendre-Gauss-Lobatto  points 
and  let  D  be  the  corresponding  differentiation  matrix  (3.10).  If 

~  _  2  1 
“  Wj  ~  N(N  +  1)  L2N(Cj)  ’ 


||  •  || is  the  Frobenius  norm.  The 
-  silt,  (3.43) 


(3.45) 
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then 


Wj  Dj  /,■ 

wkDkj , 

J  ~f~  k, 

WjDjj  = 

—WjDjj, 

j  =  1,  •  • 

,,iV-l 

27^00  = 

-i  M), 

2Dnn  = 

l/wN. 

(3.46) 


Proof  This  result  can  easily  be  verified,  keeping  in  mind  that  L^(— 1)  =  (— 1)^, 
Ln{  1)  =  1-  □ 


For  the  Chebyshev  collocation  the  following  lemma  provides  a  suboptimal  solution 
of  the  linear  least  squares  problem  (3.43). 


Lemma  3.5  Let 

Cj  =  -  COS  ,  j  =  0, . . . ,  N, 

(Chebyshev-Gauss-Lobatto  collocation  points),  let  D  be  the  correspond¬ 
ing  differentiation  matrix  (3.7).  If 

Wj  =  Wj  = 

then 

Wj  Dj  k 

27^00 

2Dnn 

The  least  squares  norms  \\W D  +  DTW  —  E\\2F,  using  the  weights  defined  by  (3.47), 
for  different  numbers  of  collocation  points  are  shown  in  Figure  3.3. 

For  the  Chebyshev  collocation  points  the  linear  least  squares  problem  (3.43)  is 
solved  numerically.  The  least  squares  norms  ||bF77  +  DTW  —  E\\2F,  using  optimal 
weights,  for  different  numbers  of  collocation  points  are  shown  in  Figure  3.4.  Figure  3.5 


JL.  j  —  n  N 

2  N  J 

2L  j  —  l  N  -  1 
TV  J  ir  •  •  -1-? 


(3.47) 


-wkDkj,  j  ^  k, 

-l/w0, 

l/wN. 


(3.48) 
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Figure  3.3:  The  least  squares  norms  ||WTD  +  DJW  —  E\\2F  for  w  given 
by  (3.47)  for  different  numbers  of  collocation  points 


Figure  3.4:  The  least  squares  norms  ||VF.D  +  DTW  —  E\\2F  for  different 
numbers  of  collocation  points 
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Figure  3.5:  min  Wi  and  max  for  different  numbers  of  collocation  points 


shows  that  the  optimal  weights  Wj,  j  =  0 , ,N,  are  positive.  The  element  wise  error 
| WD  +  DTW  —  E |  for  N  =  64  is  displayed  in  Figure  3.6. 

If  we  define 

6jk  wkDkj  -j-  WjDjk  Ejk ,  (3.49) 


j,  k  =  0, . . . ,  N,  and  use  the  identities  (3.40)  in  (3.13a),  we  obtain 


Wo 

fy(V0,  uo)T  +  J2k=oDok~kk 

+  t^  +^0 

Wj 

+  Efc=o 

WN  fy{yN,U]y)  +  Efc=0  Djk^k 

“ife  +  my(yN)  +  by(yN)T^N 


\r~^N 

2_^k=0  e0k 

2^k= 0  ejk 


y fc 
wk  ’ 

Wk  ’ 


E 


N 

k= 0 


_  Xu 
eNk^- 


j  —  1, . . . ,  N  —  1,  (3.50) 
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Figure  3.6:  Element  Wise  Error  log10  \  {WD  +  DTW  —  E)ij\  for  N  =  64. 


w0 


Wi 


We  use  the  adjoint  equations  (2.6b)  evaluated  at  cv  j  =  0, ...,1V,  and  the 
transversality  conditions  (2.6c),  to  obtain 

fy(y{co),u{c0))Tp(c0)  +  J2k= o  A)  kP(ck) 

+p(co )  +  go  =  wo  [£PN(p)(c0)  -  ftp(c0)]  , 

fy(y(cj),  m(cj))tp(cj)  +  Ef=o  Djkp(ck )  =  Wj  [ftPN{p)(cj)  -  £p(cj)]  , 

j  =  1, 

Av  fy(y(cN),u(cN))Tp(cN )  +  Ef=0  Djkp(ck ) 

-P(cjv)  +  my(y(cN))  +  by(y(cN))Tqf  =  Ar  [^-Pjv(p)(cjv)  -  |p(civ)]  • 

(3.51) 

Recall  the  definition  (3.20)  of  F^. 

Subtracting  (3.50)  from  (3.51)  leads  to  the  following  result. 


Lemma  3.6  (Consistency)  Let  p  satisfy  the  adjoint  equation  (2.6b),  let 
q0  and  qj  satisfy  the  transversality  conditions  (2.6c),  and  let  Xj,  j  = 
0,...,1V,  po,pN  satisfy  (3.19a).  If  fy{y{cj),  u{cj))  =  fy(yj,Uj),  j  = 
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0,  •  •  • ,  N,  by(y(cN ))  =  by(yN)  and  my(y(cN ))  =  my(yN),  then 
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W0 


fy(y(c0),u(c0))T  \p(cq)  -  ^  ]  +^TDok  lp(ck)  -  ^ 
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eofc  — , 

z  n 
fc= 0 


Wi 


=  Ws 
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fy{y{cj)i  u(Cj))T  (  p(cj)  -  ^  )  +  X)  DJk  (  _  ^ 

aV  A:=o  V  fc 


jtPM(Cj)  ~  ip(Ci) 


AT 


Az 


t^o  Wk 


WN 
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fy(y(cN),  U(CN))T  (  p(cN )  -  ^  Dik  f  ~  ^ 

WNJ  t^o  \  Wk 


A 


p{cN)  -  ^  J  +  my(y(cN ))  +  by(y(cN))T(qf  -  pN) 
Wn 


=  U>N 


jtPN(p)(cN)  -  jtP(cN) 


N 

E 

k= 0 


t-Nk- 


A  k 

Wk 


(3.52) 


In  the  case  of  Legendre- Gauss-Lobatto  points,  Cj,  j  =  0 , ,N,  the  weights 

~  _  2  1 
Wj  Wj  ~  N(N  +  1)  L2N(cj)  ’ 

are  optimal  and  lead  to 

E'fc  =  0,  j,k  —  0, . . .  ,N. 

In  this  case  the  adjoint  estimates  (3.18)  and  (3.39)  are  identical.  Furthermore,  in  this 
case  the  consistency  results  in  Lemma  3.3  and  in  Lemma  3.6  are  identical.  However, 
for  the  Chebyshev- Gauss-Lobatto  points,  Remark  3.8  shows  that 

\€jj\  >  0,  j  =  0,...,N, 
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and  the  numerical  results  displayed  in  Figure  3.4  indicate  that 

N  N 

^ e*  ~  °-65  > 0 

j= 0  k= 0 

Hence,  in  the  case  of  Chebyshev-Gauss-Lobatto  collocation,  the  adjoint  estimate 
(3.39)  suggested  by  the  weighting  matrix  approach  is  not  useful,  unlike  the  adjoint 
estimate  (3.18)  derived  earlier. 


3.4  Discretization  Error  for  the  Optimal  Control 


With  the  adjoint  estimation  procedure  in  place,  it  is  now  possible  to  quantify  the 
error  between  the  state 

N 


control 


yN(t)  =  ^y^t), 

i= 0 

N 


and  adjoint 


uN(t)  =  5>^(t), 

i=0 

N 


A  N(t)  =  A. Mt), 

i= 0 

computed  as  the  optimal  solution  of  the  discretized  optimal  control  problem  (3.5)  and 
the  solution  y,  u,  and  p  of  the  infinite  dimensional  optimal  control  problem  (2.5). 
Recall  that  yv  Uj,  X j,  j  =  0, N  and  p0?  Ahv  satisfy  the  weighted  discrete  adjoint 


equations 


fy(yo,  u0)Tw0X0  —  Y2k=o  DkfiWkXk  —  —yo, 

fyillji  Uj^WjXj  -  J2k=o  Dk,jWk Xk  =  0,  j  =  1, . . . ,  N  -  1,  (3.53a) 

fyiljN,  un)tu>nXn  —  J2k= o  Dk,NWkXk  =  —  by(yN)TpN  —  rny(yN), 

the  weighted- discrete  gradient  equations 


fu(Vj ,  UjfwjXj  =  0,  j  =  0, . . . ,  N, 


(3.53b) 
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and  weighted-discretized  state  equations 

{f(yvui)  - Ef=o DJ:ky^j  =  o,  j  =  o,...,N, 

yo-yo  =  0,  (3.53c) 

b(yN)  =  0. 


If  /  is  affine  linear, 


f(y(t),u(t ))  =  Fy(t)y{t)  +  Fu(t)u(t)  +  fait), 

if  m  is  quadratic, 

and  if  b  is  affine  linear 

b(y(tf ))  =  By(tf)  +  ba, 

then  the  optimality  conditions  (3.53)  can  be  written  as 


KjvXtv  =  hN, 


(3.54) 


where 


_  /  1  T  1  T  \T  \T  T  T\T 

xAf  ~  \y  0  >  •  •  ■  j  Vni  u0i  ■  ■  ■  i  UNi  )  •  •  •  )  Vo  j  Vn) 


(3.55) 


Lemma  3.7  If  y,u  satisfy  the  state  equations  (2.5b)-(2.5d),  then 


N 


k= 0 


wj[f(y{cj),uicj))  ~^2Djky{ck)]  =  raj(N,y,  1),  j=0,...,N, 

) 

y(co)  -  Vo  =  0, 

b(y(cN ))  =  0>  (3.56) 


where 


rsAN,p,g)  = 


Wj 


d 


d 


-PNip)iCj)  ~  ~ p(Cj ) 


dt 


(3.57) 


9{Cj)  LL 

Proof  This  result  follows  from  evaluating  (2.5b)  at  the  collocation  points  and  using 
the  definition  (3.20)  of  the  interpolating  polynomial.  □ 
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The  following  error  results  are  shown  for  linear-quadratic  optimal  control  prob¬ 
lems.  More  analysis  is  needed  to  extend  these  results  to  nonlinear  OCPs,  however 
that  exceeds  the  scope  of  this  thesis. 

The  following  lemma  provides  a  consistency  result  which  will  be  used  to  derive  an 
error  estimate  for  the  optimal  control  for  linear-quadratic  OCPs. 

Lemma  3.8  (Consistency  for  Linear- Quadratic  OCPs)  Let  /  be  affine 
linear, 

f{y{t),u{t ))  =  Fy[t)y(t )  +  Fu(t)u(t)  +  fa(t), 
let  m  be  quadratic, 

fn{y{tf))  =  \y{tf)TMy{tf)+mJy{tf)+ma, 
and  let  b  be  affine  linear 


Ky(tf))  =  By(tf)  +  ha 


If  y,u,p,qo,qf  are  the  solution  of  (3.5)  and  corresponding  adjoint  vari¬ 
ables  and  Lagrange  multipliers,  and  if  y0, . . . ,  y^,  u0, . . . ,  u at,  A0, . . . ,  Ajv, 
Ho,/J>n  are  the  solution  of  the  discretized  optimal  control  problem  (3.5) 
and  corresponding  weighted  Lagrange  multipliers,  then 


rp  /  nT  ( P(C o) 

Fs(c„) 


p  f  P(Ck ) 


fc= 0 


\9(ck) 

+  (<?o 


Az 


Po) 


Fy{Cj)TWj 
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Aj  j  y]  Dp-jU'l; 


k= 0 


P{Ck) 

g(ck) 


Fy(CN)TWN 


f  p(cN) 

\g(cN) 


\  N 

—  Ajv  j  —  DkNWk 

'  k= 0 


/  P(Cfc) 
V^(O-) 


+5r(g/  -  /ijv)  +  M{y{cN)  -  yN) 


ro(N,p,  g)i 
r)(N,p,g ), 
j  =  1, . . . ,  iV  —  1, 


raN(N,p,g), 


(3.58a) 
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Fuic/fvJ,  -  Aj)  =  0,  j  —  0, , . . ,  N,  (3.58b) 

“\i  Fv(cj)(y(cj)  ~  Vj)  +  Fu(cj)(u(cj)  -  Uj ) 

N 

-  Djk(y(ck)  -  yk)  =  rj  (N . ,  y,  1), 

k= 0 

j  =  0, . . . ,  N, 

y(c0)  -yo  =  o, 

B(y(cN)  -  yN)  =  0,  (3.58c) 

where  r°;(N,y,g )  and  rj(N,y,g),  j  =  0 are  defined  as  in  (3.26) 
and  (3.57)  respectively. 

Proof  The  equations  (3.58a)  were  derived  in  Lemma  3.3.  The  equations  (3.58b) 
are  obtained  by  evaluating  (2.6d)  at  c3  and  subtracting  (3.53b).  The  equations  (3.58c) 
are  obtained  by  subtracting  (3.53c)  from  (3.56).  □ 

The  first  part  of  the  following  theorem  is  an  immediate  consequence  of  Lemma 
3.8.  parts  two  and  three  follows  from  Corollary  3.1. 


Theorem  3.9  (Error  for  Linear- Quadratic  OCPs)  i.  Let  the  assump¬ 
tions  of  Lemma  3.8  be  valid.  If  is  defined  as  in  (3.55)  and  if 


x  =  y(c0)T, . . . ,  y(cN)T,  u(c0)T, u(cNy 


•  p(c0)T  p(cN)T 

’  g(co)  ’ " '  ’  g(cN)  ’ 


T  T 

%  ,Qf 


then 

||xjv  x|| 2  <  ||K^1||2||r(iV,|/,p,^)||2, 


(3.59) 
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where  Kjy  is  the  system  matrix  in  (3.53),  (3.54)  and 


/ 


r(N,y,p,g )  = 


ro(N,p,g) 

raN(N,p,g) 

0 

0 

rs0(N,y,l) 

r8N{N,y,l) 

0 

0 


\ 


(3.60) 


ii.  Let  Co, . . .  ,cjv  be  the  Chebyshev-Gauss-Lobatto  collocation  points  de- 
hned  in  Example  3.2  and  g(t )  =  l/\/l  —  f2.  If  y  and  p  are  s- times  con¬ 
tinuously  differentiable,  s  >  2,  and  er  >  0,  then  there  exists  a  constant  C 
independent  of  y,  p  and  N  such  that 


\\r{N,y,p,g)\\2 

<  cn 2  6  (IblUff  +  IblUa  +  Ib/fi'lUg) 

+C'Ar-ff||p||Siff  oma^  Uj/gW^g-  (3.61) 

iii.  Let  Co, . . . ,  cat  be  the  Legendre-Gauss-Lobatto  collocation  points  de¬ 
fined  in  Example  3.3  and  g(t)  =  1.  If  y  and  p  are  s-times  continuously 
differentiable,  s  >  5/2,  then  there  exists  a  constant  C  independent  of  y, 
p  and  N  such  that 

\\r(N,y,p,  1)||2  <  CN^2-9{\\y\\a>1  +  \\p\\atl).  (3.62) 


To  obtain  an  error  estimate,  one  needs  a  stability  result  that  guarantees  the  uni¬ 
form  boundedness  of  1 1 K^1 1 1 2 -  Such  a  result  is  not  yet  known. 
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Example  3.10  Consider  Example  2.13.  Applying  a  Legendre  pseu- 
dospectral  discretization  to  this  problem  yields  the  results  shown  in  Figure 
3.7.  The  numerical  results  indicate  that  the  solutions  of  the  discretized 
problem  converge  quickly  to  the  solution  of  the  infinite  dimensional  prob¬ 
lem.  The  lower  right  plot  in  Figure  3.7  also  shows  that  1 1 K ^  1 1 2  increased 
significantly  as  N  increases. 


Figure  3.7:  Error  vs.  N  for  Linear-Quadratic  Optimal  Control  Problem 
in  Mayer  Form  Using  Legendre  Pseudospectral  Collocation- 
Top  Left:  i 2  State  Error  -  Top  Right:  i 2  Control  Error 
-  Bottom  Left:  i2  Adjoint  Divided  by  Weighting  Function 
Error  -  Bottom  Right:  Norm  of  System  Matrix  Inverse 
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Example  3.11  Again  consider  Example  2.13.  Now  the  Chebychev  pseu- 
dospectral  discretization  is  applied  to  this  problem.  The  numerical  results 
are  shown  in  Figure  3.8.  The  error  between  the  solutions  of  the  discretized 
problem  and  the  solution  of  the  infinite  dimensional  problem  decays  much 
slower  than  in  Example  3.10.  Especially  the  error  XN  —  p/g  for  given  N  is 
much  larger  than  in  Example  3.10.  We  also  observe  that  for  the  Chebychev 
pseudospectral  discretization  || K^1 1| 2  is  larger  and  increases  more  rapidly 
as  N  increases  than  in  Example  3.10. 


Figure  3.8:  Error  vs.  N  for  Linear-Quadratic  Optimal  Control 
Problem  in  Mayer  Form  Using  Chebyshev  Pseudospectral 
Collocation-  Top  Left:  £ 2  State  Error  -  Top  Right:  £2  Control 
Error  -  Bottom  Left:  £2  Adjoint  Divided  by  Weighting 
Function  Error  -  Bottom  Right:  Norm  of  System  Matrix 
Inverse 
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From  the  numerical  results  in  Examples  3.10  3.11  it  is  questionable  whether  one 
can  prove  that  || K^r1 1| 2  is  bounded.  The  numerical  results,  however,  indicate  that 
even  if  || K^1 1| 2  is  not  bounded  it  grows  slower  than  ||r(!V,  y, p,  g)\\2  decreases.  In  such 
a  case,  convergence  of  the  solutions  to  the  discretized  problems  can  be  guaranteed, 
but  the  rate  of  convergence  is  less  than  one  would  expect  based  on  the  consistency 
results  alone.  It  is  also  not  known  whether  and,  if  so  how,  the  growth  in  HK^r1  H2  is 
related  to  the  increasing  condition  number  of  the  constraint  Jacobians  reported  on  in 
Example  3.1. 

3.5  Numerical  Equivalence  of  Bolza  and  Mayer  Forms 

In  this  section  the  numerical  difference  between  the  Bolza  form  OCP  (2.1)  and  the 
Mayer  form  OCP  (2.9)  is  addressed.  In  [17]  the  argument  is  made  that  for  the 
Legendre  pseudospectral  method  the  quadrature  rule  used  to  compute  the  integral 
in  (2.1a)  is  equivalent  to  the  resulting  auxiliary  discrete  adjoint  equations  in  the 
transformed  problem  (2.9).  This  would  imply  that  for  the  Legendre  pseudospectral 
method,  a  direct  transcription  of  either  problem  leads  to  the  same  numerical  solution. 
It  will  be  shown  that  this  assertion  is  not  quite  correct  and  that  solving  the  discretized 
OCP  in  Mayer  and  Bolza  forms,  respectively,  yield  results  that  merely  converge  to 
the  same  solution  as  IV  — >  00.  Recall  the  Bolza  form  optimal  control  problem 

min  m{y{  1))  +  J  £(y(t),u(t))dt,  (3.63a) 

=  f(y(t),u(t)), 

y(- 1)  =  yo, 

b(y(  1))  =  0. 


s.t. 


(3.63b) 

(3.63c) 

(3.63d) 
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Recall  the  transformed  Mayer  form  OCP 


min 

s.t. 


™(i/(i))  +  *(!), 

(3.64a) 

=  f(y(t),u(t)), 

(3.64b) 

Jtz(t)  =  £(y(t),u(t))  | 

(3.64c) 

J/(-l)  =  J/o, 

(3.64d) 

^(— 1)  —  o, 

(3.64e) 

%(!))=  0. 

(3.64f) 

To  compare  the  discrete  solutions  to  (3.63a)  and  (3.64)  it  is  necessary  to  look 
at  the  discrete  optimality  systems  of  each.  The  weighted  discrete  Lagrangian  from 
(3.16)  for  (3.64)  can  be  written  as 

Lw(y,  z,  u,  A,  7,  /i0,  nN,  z/0)  =  m(yN)  +  b(yN)TyN  +  (y0  -  y0)T  y0  +  z%v0 


N  N 

+Wjlj  [E  £{Vj >  ^  ^  ^ j,kZk\ 


3=0 

N 


k= 0 
N 


+wj\j 


fiy*’ u j)  -  D3,kVk]  •  (3-65) 


j=0  k= 0 

Differentiating  (3.65)  with  respect  to  the  yj s  and  setting  it  equal  to  zero  yields 
fy(y o,  u0)Tw0X0  +  £y(yo,  u0)w0rro  —  J2k=o  Dk,owkXk  =  —  yo, 

fy(yji  UjfwjXj  +  £y(yj,  Uj)Wjlj  -  Y*=o  Dk,jWkXk  =  0,  j  =  1, . . . ,  N  —  1, 

fy (jJn j  un)twnXn  +  £y{yNi  un)wn^n  —  Ylk= o  Dk,NWkXk  =  —by(yN)TyN  —  rny(yN) 

(3.66a) 

Differentiating  (3.65)  with  respect  to  the  z/s  and  setting  it  equal  to  zero  yields 


SfcLo  DkA^klk  —  —  Dn 

Ei=o  Dkjwk 7fc  =  0,  j  =  1, . . . ,  TV. 


(3.66b) 
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Differentiating  the  Lagrangian  (3.65)  with  respect  to  the  Uj’s  and  setting  the  deriva¬ 
tives  to  zero  gives 


fu(Vj,  UjfwjXj  +  =  0,  j  =  0, . . . ,  N. 


(3.66c) 


Differentiating  the  Lagrangian  (3.65)  with  respect  to  the  A/s,  po  and  /qv  and  setting 
the  derivatives  equal  to  zero 


Wj 


=  0,  j  =  0, . . . ,  N, 


f  ( Uj  j  uj)  Y2k=< o  Dj,kyk 

yo-yo  =  0,  (3.66d) 

Kvn)  =  0. 

Differentiating  the  Lagrangian  (3.65)  with  respect  to  the  A/s,  /iq  and  /qv  and  setting 


the  derivatives  equal  to  zero  as  well 


ID, 


i  (yj ,Uj)  'Yhk= 0  D j,kzk 


=  0,  j  =  0, . . . ,  N, 


(3.66e) 


=  0. 


Alternatively,  the  OCP  (3.63)  can  be  solved  directly  in  Bolza  form.  In  this  case 
the  integral  term  is  approximated  by 

[  £(y(t),u(t))dt  =  I  ^§-£(y(t),u{t))dt 


>- i 


l-i  9(t) 


N 


Using  (3.67),  the  weighted  discrete  Lagrangian  for  (3.63)  can  be  written  as 


(3.67) 


N 


Lw(  y,u,A,/x0,/qv)  =  m(yN) +  ^-1-i(yj,uj) 

j=o  d\cj> 

+b(yN)T/iN  +  (y0  -  y0)Tdo 

N  N 

+wjXJ  [  /(%■»  ui)  -  Di*yk]  ■ 


j= 0 


k= 0 


(3.68) 
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Differentiating  (3.68)  with  respect  to  the  y/'s  and  setting  it  equal  to  zero  yields 

fv(y o,  u0)Tw0X0  +  £y(y0 ,  «o)^y  -  £jJL o  DkfiwkXk  =  -y0, 

fy(lJjiuj)  wj^j  +  ^yilJji  ui)~g(cj)  ~  Sfc= o  Dk,jWk\k  =  0,  j  =  1, . . . ,  N  —  1, 

fy(]jN,  un)twnXn  +  £y(yN,  ~  J2k= o  D k,NwkXk  =  —by(yN)T yN  —  my(yN). 

(3.69a) 

Differentiating  the  Lagrangian  (3.68)  with  respect  to  the  u/ s  and  setting  the  deriva¬ 
tives  to  zero  gives 


U)  ■ 

fu(Vj,  u^WjXj  +£u(yj,uj)-J-  =  0,  j  =  0, . . . ,  N.  (3.69b) 

9\cj) 

Differentiating  the  Lagrangian  (3.68)  with  respect  to  the  A/s,  /i0  and  y^  and  setting 
the  derivatives  to  zero  gives 


Wn 


/('//•  ".i 


2-^/k=\ 


o  Djtkyk 


Vo  ~  Vo 


0,  j  =  0, . . . ,  N, 

0, 


(3.69c) 


Kvn)  =  0. 

Lemma  3.9  Let  yM,uM,  XM ,  ,  yff ,  y^  and  u^'1  be  solutions  to  the 

weighted  discrete  optimality  system  (3.65)  corresponding  to  the  trans¬ 
formed  Mayer  form  OCP  (3.64).  Let  yB,uB,XB,yB  and  yB  be  solutions 
to  the  weighted  discrete  optimality  system  (3.68)  corresponding  to  the 
Bolza  form  OCP  (3.63).  We  define 


XM  - 

Xjv  — 

(to; T)T.. 

■  ■  ,{Vn)T  Auv)T ,  ■ 

■  ■ ,  (4T,  (A0M)T,  •  •  • ,  (A]v)T,  (y^f,  (/4T) 

and 

xB  - 

(W)T... 

■  • ,  (^f,  («?)T,  ■  ■  ■ 

* (4)r.  (A^)T, •  •  • ,  (A^)T,  ( yo)T ,  ( /4)T )  , 

to  be  the  numerical  solutions  to  (3.64)  and  (3.63)  respectively.  If  the 
weighted  discrete  optimality  systems  (3.65)  and  (3.68)  are  sufficiently  sta¬ 
ble,  then  the  error  between  these  solutions,  ||x:(/  —  x^||2,  will  converge  to 
zero  as  N  — >  oo. 


53 


Proof  This  result  is  easily  verified  using  the  adjoint  estimation  Lemma  3.2  and  the 
optimality  conditions  (2.10)  in  Section  2.2,  7  — >  r(t)/g(t)  =  l/g(t),  where  r(t)  is  the 
true  auxiliary  adjoint.  □ 

It  is  evident  that  the  systems  (3.66)  and  (3.69)  are  not  equivalent.  This  is  because 
the  estimated  auxiliary  adjoint, 

N 

7  N(t)  = 

k= 0 

will  only  converge  to  its  true  solution  1  / g(t)  as  N  — >  00.  In  order  for  these  systems  to 
be  equivalent,  the  auxiliary  adjoint  7 N  would  have  to  be  equal  to  its  true  solution,  1, 
for  all  N.  This  notion  is  reinforced  by  the  following  example  to  conclude  this  section. 

Example  3.12  Consider  the  example  problem  (2.2)  in  Mayer  from  (2.13). 

A  Legendre  pseudospectral  discretization  to  this  problem  is  applied.  Secondly, 
consider  the  example  problem  (2.2)  in  Bolza  from  (2.11).  Again  a  Legendre 
pseudospectral  discretization  to  this  problem  is  applied.  Taking  the  £2 
norm  error  between  the  state,  control,  and  adjoint  for  each  N  yields  the 
results  shown  in  Figure  3.9.  Notice  that  the  behavior  described  in  Lemma 
3.9  is  exhibited.  The  solutions  are  never  identical,  but  converge  to  the 
true  solution  a,s  N  —>  00. 
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Figure  3.9:  Q  Error  between  Mayer  Form  Problem  and  Bolza  Form 
Problem  vs.  N  for  Linear-Quadratic  Optimal  Control 
Problem  Using  Legendre  Pseudospectral  Collocation-  Top 
Left:  State  Error  -  Top  Right:  Control  Error  -  Bottom 
Middle:  Adjoint  Divided  by  Weighting  Function  Error 


55 


3.6  Extension  to  Multiple  Subintervals 

Much  of  the  work  done  in  this  chapter  applied  the  pseudospectral  method  on  one 
interval  [—1,1].  Extensions  to  multiple  intervals  are  very  important  for  many  prob¬ 
lems.  Our  error  bound  in  Theorem  3.9  indicates  that  the  discretization  error  between 
computed  solution  and  true  solution  depends  on  the  smoothness  of  the  state  and 
of  the  adjoint.  The  smoothness  of  the  state  depends,  among  other  things,  on  the 
properties  of  the  right  hand  side  function  /  in  the  governing  dynamics.  For  problems 
with  piecewise  continuous  right  hand  sides  (e.g.,  due  to  change  of  mass  in  launch 
problems,  or  dne  to  piecewise  constant  controls),  it  is  important  to  introduce  mul¬ 
tiple  subintervals.  Another  potential  benefit  of  using  pseudospectral  methods  along 
multiple  intervals  is  to  take  advantage  of  sparsity.  Indeed  the  optimality  system  for  a 
pseudospectral  method  along  many  subintervals  will  be  very  sparse  relative  to  a  pseu¬ 
dospectral  method  applied  on  one  interval.  The  benefit  is  that  proper  exploitation  of 
sparsity  may  improve  solution  time. 

The  pseudospectral  method  can  easily  be  extended  to  multiple  subintervals.  To 
accomplish  this,  the  collocation  points 

C0  =  — 1,  Ci,  ...  ,  Cat-i  G  (  —  1, 1),  Cat  =  1 

are  again  used.  At  this  point  it  is  more  useful  to  consider  the  OCP  (2.5)  on  the 
interval  [to,tf]. 

Remark  3.13  The  time  interval  shift  can  be  accomplished  by  the  fol¬ 
lowing  identity.  Let  t(c)  G  [tf,t0\  be  the  mapping 

t(c)  =  ((tf  ~  to)c  +  tf  +  t0)/2. 

By  the  chain  rule,  we  have  that 

^2/(*(c))  =  ^Wc))^(c)  =  tf  2  tQ/(^(c))X*(c)))- 
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However,  for  the  remainder  of  this  section,  t  will  not  be  written  as  an  explicit  function 
of  c  and  the  more  convenient  notation 

jty{t)  =  tf  2  to f(y(t),u(t)), 

will  be  used. 

The  interval  [t0,  tf]  is  subdivided  into  /  subintervals  [ti,  ti+ 1],  i  —  0, —  1,  with 


t0  <ti  <  ...  <t!  —  tf. 


We  define  hi  =  U+i  — ti .  The  state  y  is  approximated  by  a  piecewise  polynomial  yh,N . 
The  restriction  of  yh,N  onto  [ti,  ti+ 1],  i  —  0, —  1,  is  denoted  by  y-  'N  and  written 
as 


h.N 


N 

tt)  =  Y.y 

j= 0 


Vb 


-1  +  2 


t  -  ti 
hi 


(3.70) 


The  collocation  discretization  of  the  optimal  control  problem  (2.5)  is  given  by 


min  m(y/_iiJv), 
s.t. 


t  y,o  N  ( 


D 


hi. 

2 


\  Vi,N  / 


f(yi,o,Uifi)  ^ 

V  f  (,Ui,Ni  W*,Jv)  J 


Ui.N  ~  Vi+ 1,0, 


1/0,0  ~  VOi 


(3.71a) 

i  =  0, ...,/  —  1,  (3.71b) 

i  =  0, ...,/  —  2,  (3.71c) 

(3.71d) 


b(yi~i,N)  —  0. 


(3.71e) 
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The  weighted  Lagrangian  corresponding  to  (3.71)  is  given  by 

Lw(y,u,X,jl,no,liN)  =  m(yI-ijN)  (3.72) 

1-1  N  ^  N 

Dj^kVi,k  J 

i= 0  j= 0  k= 0 

1-2 

+  2^  ~  Vi+ i,o] 

i= 0 

+b{yi-i,N)T UN  +  (2/0,0  -  j7o)V  (3.73) 

Differentiating  the  Lagrangian  (3.72)  with  respect  to  the  y^' s  and  setting  the  deriva¬ 
tives  to  zero  gives  the  weighted-discrete  adjoint  equations  on  multiple  intervals 

^fy(y 0,0)  uO,o)TwO^O,0  ~  EfcLo  Dk,0wk^0 ,k  +  A*0 j  —  0? 

tffyiVoji  uo,j)TWjX0,j  ~  EfcLo  Dk,jWk Xo,k  =  0,  j  =  1, . . . ,  N  -  1, 
^fy(yo,N,  Mo,w)TW7vAo,Ar  —  EfcLo  Dk,NU>k\o,k  +  fio  =  0, 

(3.74a) 

on  the  first  subinterval  i  =  0, 

~2  fy  (2/*,0  5  ^h,o)  ^oAi,0  ^  A — n  Pk.Q^k  Aj.fc  /b— 1  0, 

~2  fy(yi,ji  ^ i,j )  ‘Wj^i,j  ^2k= 0  ^k,j^kXi,k  0,  j  1,  .  .  .  ,  iV  1, 
\fy{yi,Ni  Ui}N)TWN\iyN  —  EfcL o  Dk)NWkXi,k  +  P*i  =  0, 

(3.74b) 

for  i  =  1, ...,/  —  2,  and 

\  fy{y  1-1,0,  W/-1,o)TWoA/_1)0  —  Efc^O  DkflWk\l-l,k  —  ftl-2  =  0, 

~2  fyiVl—  l,j;  Ul—  Ij)  WjXj—ij  Efe=0  Dk,jWjXl— Itk  0,  J  1 ......  A"  1, 

\  fy(yi-l,N,  Ui-1^n)TWnXi-1^n  ~  Efc= 0  Dk,NWkXl-l,k 

+by(yI-i)N)TnN  +  my(yI-^N)  =  0, 


on  the  last  subinterval  i  =  /  —  1. 


(3.74c) 
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The  next  lemma  shows  that  the  adjoint  variable  p  divided  by  the  weighting  func¬ 
tion  g  satisfies  the  weighted-discrete  adjoint  equations  along  multiple  subintervals 
(3.19a)  with  an  error  that  is  dependent  on  the  true  adjoint  p,  the  weighting  function 
g  and  N. 


Lemma  3.10  If  p,  q0  and  qg  satisfy  the  adjoint  equation  (2.6b)  and  the 
transversality  conditions  (2.6c)  then 


0  r  xtP(^0,0 


N 


~^fy(yo,0,Uo,0 


k= 0 


P(t0,k 


9(c o)  £-7  k'°  g(ck)  +  q° 


=  rZ0(N,p,g), 


ho  r  i 

-jfy{yo,j,uo,r 
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on  the  first  subinterval  i  =  0, 
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for  i  —  1, —  2,  and 


hi  r  /  i,o) 

yM9/-i,o,«/-i,o) 


TV 


E°* 

k=0 


p{tl-l,k) 


tM  t'-k'°^r~p{t,-2'N) 


=  fj_i)0  (N,p,g), 


N 


=  ri-ij(N,p,  g),  j  =  1, 

hi  ,  ,  stP^i- i,n)  n  p(^/-i,fc) 

+&?/(l//_i)Ar)Tg/  +  my(yI-^N) 

=  rj_1>N{N,p,g), 
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on  the  last  subinterval  ?'  =  /  —  !,  where 
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(3.77) 


Proof  The  result  is  a  direct  extension  of  Lemma  3.2.  It  is  obtained  by  using 
equation  (3.23)  and  the  fact  that  p,  qo  and  qf  satisfy  the  adjoint  equation  (2.6b)  and 
the  transversality  conditions  (2.6c).  Then  p/ g,p}  q0  and  qf  are  inserted  into  (3.74)  for 
\,jl,Po  and  pn  respectively  to  obtain  (3.75).  □ 
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Chapter  4 

International  Space  Station  Momentum  Dumping 

Problem 

The  Legendre  Pseudospectral  method  described  in  the  previous  chapters  is  now  ap¬ 
plied  to  a  realistic  optimal  control  problem.  This  chapter  describes  formulation  and 
solution  of  the  International  Space  Station  momentum  dumping  problem.  One  version 
of  this  problem,  where  a  continuous  control  is  considered,  lends  itself  to  the  appli¬ 
cation  of  the  Legendre  Pseudospectral  method  on  one  interval  while  other  versions, 
where  piecewise  constant  controls  are  considered,  lend  themselves  to  the  Legendre 
Pseudospectral  method  using  multiple  subintervals.  In  each  case,  the  problem  is 
stated,  then  transcribed  into  a  nonlinear  program  and  solved  using  standard  nonlin¬ 
ear  programming  techniques.  Numerical  results  for  each  problem  scenario  are  given. 

4.1  Background 

Spacecraft  attitude  control  is  usually  provided  by  momentum  devices  such  Control 
Moment  Gyroscopes  (CMGs)  or  reaction  wheels,  as  they  do  not  require  consumables. 
However,  the  momentum  of  such  devices  is  limited  and  when  this  limit  is  reached  the 
device  is  termed  saturated.  In  this  situation,  ‘controllability’  is  lost  along  the  mo¬ 
mentum  saturation  direction.  Recovering  full  three  degree-of-freedom  control  requires 
desaturating  the  momentum  device. 

The  usual  approach  to  desaturate  accumulated  momentum  is  to  use  an  addi¬ 
tional  device.  Examples  are  mass  expulsion  devices,  magnetic  dipoles  which  interact 
with  the  Earth’s  magnetic  field,  and  rotating  solar  arrays  which  interact  with  solar 
radiation  pressure  [22], [50].  Mass  expulsion  devices  require  the  use  of  consumable 
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propellant,  which  has  finite  lifetime  and  is  expensive  to  get  to  orbit  or  replenish. 
Magnetic  dipoles  are  electromagnets,  which  generate  a  torque  on  the  spacecraft  by 
their  dipole  interaction  with  the  Earth’s  magnetic  field.  Disadvantages  are  that  the 
Earth’s  magnetic  field  is  not  well  known  and  hence  may  require  the  use  of  magne¬ 
tometers  to  measure  it,  that  it  can  be  affected  by  sun  spots  or  magnetic  storms,  and 
that  it  varies  with  orbit  location  thus  restricting  the  amount  and  direction  in  which 
momentum  can  be  unloaded.  Roll  and  pitch  momentum  is  unloaded  near  the  mag¬ 
netic  poles,  while  roll  and  yaw  momentum  is  unloaded  near  the  geomagnetic  equator 
[32], [51].  Further,  the  use  of  dipoles  generates  an  additional  magnetic  field  on  the 
spacecraft,  which  may  affect  other  sensors  or  devices.  Solar  pressure  based  methods 
require  the  use  of  modulating  surfaces  such  as  solar  arrays.  This  sacrifices  electrical 
power,  creates  mechanical  lifetime  issues  clue  to  wear  and  tear  and  increases  the  risk 
of  drive  failure. 

An  alternative  is  to  use  the  momentum  devices  to  appropriately  maneuver  the 
spacecraft  in  a  disturbance  field  such  that  accumulated  momentum  can  be  removed 
[51].  Since  most  environmental  disturbances  on  the  spacecraft  are  a  function  of  its 
attitude,  the  accumulated  momentum  due  to  navigating  in  such  a  disturbance  field 
is  path  dependent.  Performing  an  attitude  maneuver  over  a  pre-selected  trajectory 
can  result  in  a  lower  final  momentum  state  than  one  with  which  the  vehicle  started. 
The  advantage  of  this  approach  is  that  it  does  not  require  any  additional  hardware 
or  specialized  software.  Hence,  it  can  be  applied  to  any  existing  vehicles  that  use 
momentum  devices.  In  general,  the  method  provides  momentum  unloading  in  all 
axes  and  does  not  require  preferred  orbit  locations.  Gravity  gradient  and  to  a  lesser 
extent  aerodynamic  torques  are  well  defined  and  better  known  than  Earth’s  magnetic 
field.  This  approach  can  also  be  used  as  a  contingency  operational  mode  for  spacecraft 
which  use  other  actuators  for  momentum  dumping  purposes,  consequently  increasing 
the  operational  lifetime  of  satellites  already  in  orbit. 


62 


4.2  Rotational  Dynamics 

The  equations  of  motion  described  in  this  section  can  be  found  in  [32],  The  attitude 
dynamics  of  a  rigid  body  in  a  circular  orbit  are  given  as 

J—u(t)  =  Td  —  oj[t)x  ( Ju(t)  +  h{t))  —  u(t),  (4.1a) 

where  lv  :  M  i— >•  M3  is  the  angular  velocity  of  the  spacecraft  with  respect  to  an  inertial 
reference  frame  measured  in  rad/sec.  The  remaining  terms  are  h  :  M  i— >  M3  the 
angular  momentum  of  the  Control  Moment  Gyroscopes  (CMGs)  measured  in  ft-lbs- 
sec,  «  :  R  w  l3  the  control  torque  measured  in  ft-lbs,  Td  G  M3  the  external  disturbance 
torque,  J  G  M3x3  the  inertia  matrix  of  the  spacecraft  measured  in  slugs-ft2.  All  terms 
are  evaluated  with  respect  to  the  spacecrafts  fixed  body  reference  frame,  see  Figure 
4.1.  The  skew-symmetric  cross  product  operator  x  is  is  given  as 

O3  Cl  2  ^ 

0  —  Oi 

a  1  0  j 

External  disturbance  torques  may  take  many  forms.  Examples  are  gravity  gradient 
torque,  aerodynamic  torque  and  magnetic  torque.  This  model  considers  only  gravity 
gradient  torque  Tgg  which  is  given  by 

Tgg  =  3uj2orbC^  J C3,  (4. lb) 

where  uorb  =  0.0011  rad/sec  is  the  orbital  for  the  current  altitude  and  C3  is  the 
third  column  of  the  rotation  matrix  which  rotates  any  vector  in  the  local  vertical 
local  horizontal  (LVLH)  reference  frame  into  the  spacecrafts  body  reference  frame. 
It  is  assumed  that  all  other  external  torques  are  small  relative  to  those  modeled  and 
therefore  negligible. 

The  attitude  kinematics,  using  a  quaternion  formulation,  are  given  as 


1  0 
0.3 

V  -«2 


-  Wo). 


(4.1c) 
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where  q  :  M  i— >  M4  is  the  attitude  quaternion,  T  :  M4  i— »•  M4x3  is  given  by 

(  -92 (t)  - q3(t )  -g4(f) 

N  1  9i(f)  “94 (*)  93(0 

T(<?)  =  0 

94  (t)  9i(*)  -92  (t) 

^  -<?s(t)  92 (t)  Qi(t) 

and  uj0  G  M3  is  the  constant  orbital  rate  for  the  LVLH  reference  frame,  assuming  a 
circular  orbit.  The  u>(t)  —u0  term  represents  the  relative  angular  rate  with  respect  to 
the  LVLH  reference  frame,  therefore  the  quaternion  q  computed  from  (4.1c)  describes 
the  attitude  of  the  spacecraft  with  respect  to  the  LVLH  reference  frame.  It  is  standard 
for  the  control  variable  to  enter  into  the  dynamics  through  a  control  law.  For  now  it 
is  assumed  that  the  CMGs  are  controllable  directly,  resulting  in  the  control  law 

jtKt)  =  u(t).  (4. Id) 

Figure  (4.1)  shows  the  aforementioned  reference  frames  as  they  relate  to  the  Earth, 
the  LVLH  orbit  and  the  space  station. 

Additional  consideration  must  be  given  to  the  attitude  kinematics  equation  (4.1c) 
because  a  quaternion  must  always  have  a  unit  norm.  Therefore  a  path  equality 
constraint  must  be  added 

||9(t)||2  =  l,  Vie  [t0,*/].  (4.2) 

This  is  typically  a  difficult  constraint  to  satisfy.  For  simulations,  the  standard  pro¬ 
cedure  in  simulations  is  to  divide  the  current  quaternion  by  its  magnitude  during 
each  step  of  numerical  integration.  In  this  thesis  we  use  Legendre  Pseudospectral 
collocation  to  discretize  the  dynamics  in  the  optimal  control  problem.  Since  equality 
constraints  are  only  enforced  at  the  collocation  points,  it  should  not  be  expected  that 
the  unity  norm  constraint  will  be  satisfied  in  the  infinite  dimensional  sense.  Such  a 
constraint  violation  results  in  a  solution  which  has  no  physical  meaning.  Therefore  it 
is  beneficial  to  use  Euler-Rodriguez  parameters  which  are  defined  to  be  [32] 

r(t)  =  (92 (t)  9s (t)  1s(t))T- 


(4.3) 
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Figure  4.1:  Earth’s  Inertial  Reference  Frame,  Local  Vertical  Local 
Horizontal  Reference  Frame,  International  Space  Station’s 
Body  Reference  Frame  (ISS  Assembly  12 A  shown) 


Note  that  (4.3)  is  not  defined  when  qi(t)  =  0,  which  is  equivalent  to  a  180°  rotation. 
This  corresponds  to  attitudes  which  are  assumed  to  not  occur  along  the  optimal 
trajectory.  Using  the  representation  (4.3)  does  not  require  the  path  equality  constraint 
(4.2).  Via  Euler- Rodriguez  parameters,  (4.1c)  can  be  converted  to 

=  \(r(t)r(t)T  +  1  +  r(t)x)(u(t)  -  ua).  (4.4) 


The  rotation  matrix  in  (4.1b)  can  be  computed  as 

2 


C  =  I  + 


rxrx  —  rx' 


1  +  rTr 

The  resulting  attitude  dynamics  for  the  space  station  are  given  as 


=  rgg(r)  -u(t)x(Ju(t)  +  h(t))  -u(t) 
ftr(t )  =  \{r(t)r(t)T  +  I  +  r{t)x)(u(t)  -  ua(r)) 


&Kt)  =  u{t). 


(4.5) 
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4.2.1  International  Space  Station  Assembly  Stage  12 A 

This  thesis  considers  International  Space  Station  stage  12A,  which  was  originally 
scheduled  for  launch  in  December  2002.  This  assembly  makes  the  following  additions 
[2]: 

•  Delivers  second  port  truss  segment  (P3/P4  truss)  to  attach  to  first  port  truss 
segment  (PI  truss). 

•  Central  cooling  radiators,  delivered  earlier  on  flights  9A  and  11  A,  are  deployed 
from  first  starboard  (SI  truss)  port  (PI)  truss  segments. 

•  Exterior  attachments  for  Brazilian  Unpressurized  Logistics  Carriers  (ULCs)  are 
delivered. 

The  inertia  matrix  J  for  space  station  assembly  stage  12A,  shown  in  Figure  4.2,  is 
given  in  Table  4.1  [40]. 


2.8070  x  107  4.8225  x  105  -1.7168  x  107 

4.8225  x  105  9.5145  x  107  6.0260  x  104 

-1.7168  x  107  6.0260  x  104  7.6594  x  107 

Table  4.1:  ISS  12A  Inertia  Matrix  [slugs- ft2] 

Due  to  physical  limitations,  the  CMGs  must  not  reach  a  certain  momentum  mag¬ 
nitude  threshold  because  they  will  become  saturated.  This  saturation  limit  can  be 
found  in  [39]  to  be  10000  ft-lbs-sec.  This  leads  to  the  path  inequality  constraint 

\\h(t)\\2  <  hmax,  Vte[t0,tf]%  (4.6) 

where  hmax  =  10000. 
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Figure  4.2:  International  Space  Station  Assembly  12A 
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4.2.2  Boundary  Conditions 

The  boundary  conditions  for  (4.5)  are  chosen  such  that  the  spacecraft  is  at  a  Principal 
Axis  (PA)  attitude  initially  and  travels  to  a  Torque  Equilibrium  Attitude  (TEA).  A 
PA  attitude  is  an  attitude  for  which  the  gravity  gradient  torque  is  zero.  A  PA  is  a 
common  rest  attitude  for  a  spacecraft  such  as  the  space  station.  The  PA  attitude 
associated  with  the  inertia  matrix,  in  Table  4.1,  is  shown  in  Table  4.2. 


v(to) 

-9.54  x  10~6 

-1.14  x  10“3 

5.35  x  10“6 

r(to) 

3.00  x  10"3 

1.53  x  HT1 

3.83  x  10“3 

Table  4.2:  Attitude  and  Rate  Corresponding  to  Principal  Axis 

A  TEA  is  a  special  attitude  for  which  the  right-hand  side  of  the  differential  equa¬ 
tions  (4.5)  are  all  zero  when  no  control  is  exerted  in  the  vehicle.  Finding  and  reaching 
a  TEA  reduces  to  a  final  time  boundary  condition  where  a  root-finding  problem  to 
find  u>,  r  and  h  such  that  the  right  hand  side  of  the  differential  equation  is  equal  to 
zero  when  u  =  0.  A  TEA  corresponds  to  attitudes  that  can  be  held  indefinitely.  This 
is  a  desirable  attitude  because  when  a  spacecraft  is  at  a  TEA,  it  does  not  require 
attitude  control  devices  to  stay  at  that  attitude.  In  other  words,  the  final  state  of  the 
system  can  be  maintained  indefinitely  with  zero  control  effort. 

Additionally,  an  initial  value  for  angular  momentum  must  be  specified.  This  value 
can  change  from  one  simulation  scenario  to  the  next,  so  this  condition  is  somewhat 
less  strict  as  the  PA  and  TEA  requirements.  For  problems  considered  in  this  thesis, 
the  initial  value 

h(t0)  =  (5000,  5000,  5000)T  ft-lbs-sec 


was  used. 


These  boundary  conditions  can  be  written  compactly  as 


u(t0)  —  Uq, 

r(t0)  =  To, 
h(t0 )  =  h0, 
and 


b{w{tf),r{tf),h{tf ))  = 

(  J~\Tgg{tf)  -uj{tf)x{Juj(tf)  +  h(tf))) 

V  §(r(*/M*/)T  +  1  +  r(tf)*)(w(tf)  -  u0(r)) 


where  u;0,  r0  and  h0  are  given  by  the  above  initial  conditions. 
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4.3  ISS  Momentum  Dumping  Problem  with  Continuous  Control 

From  a  modeling  perspective,  the  simplest  version  of  the  space  station  momentum 
dumping  problem  is  posed  such  that  the  control  in  (4.5)  is  a  continuous  function  on 
the  interval  [to,tf].  The  optimal  control  problem  can  be  stated  as 

min  \\h(tf)\\2 
s.t. 


=  Tgg(r)  -  LV(t)X(ju(t)  +h(t))  ~u(t),  te[t0,tf], 

fAt) 

=  l(r(t)r(t)T  +  I  +  r(t)x)(u(t) -u0(r)),  t  e  [t0,tf}, 

ih(t) 

=  u(t ),  te[t0,tf], 

\m\\2 

—  bimax:  t  G  [tojtj], 

u(to) 

=  k>o, 

r(to) 

=  h0, 

h(t0) 

=  h0, 

b(w(tf),r(tf),h(tf )) 

=  o, 

(4.8) 


where  b,  Cdo,  fo  and  h0  are  given  by  (4.7).  This  problem  is  posed  on  the  interval 
[to,tf]  =  [0,1800]  sec  [38].  The  initial  data  uq,  L0  for  the  attitude  and  the  angular 
rate  were  chosen  to  be  the  principal  axis  from  Table  4.2.  The  initial  value  for  the 
angular  momentum  was  chosen  to  be  h0  =  (5000,  5000,  5000)T,  where  1 1 7z0 1 1 2  is 
close  to  hmax  to  make  a  desaturation  maneuver  meaningful.  The  final  time  boundary 
condition  is  defined  in  (4.7). 

Since  the  control  is  continuous  on  the  interval  [f0,£/]  the  Legendre  Pseudospectral 
method  with  one  time  interval  can  be  applied  to  discretize  this  problem.  Applying 
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this  direct  transcription  to  (4.8)  results  in  the  following  NLP 


min  \\hN\\2 

s.t. 

N 

tf-toJ^DjkUk 

J  k= 0 

=  Tgg(rj ) 

— 

N 

tf~toPoDjkrk 

=  2  (rn 

+ 

N 

tf~toPoDjkhk 

=  j 

= 

WMU 

L  hmax , 

j 

Uq 

—  ^0, 

To 

=  b0, 

ho 

=  ho, 

b(tL>N,  T TV,  hjv) 

=  o, 

’j  {Juj  +  hj)  -  Uj,  j  =  Q,...,N, 

+  rf)(uj  -  oj0(rj)),  j  =  0,...,N, 

=  0 


(4.9) 


where  the  optimization  variables  are  Uj,  rj,  hj ,  and  Uj,  j  —  0, . . . ,  N.  Note  that  the 
path  inequality  constraint  is  only  enforced  at  the  collocation  points,  ft  is  assumed  that 
doing  so  will  result  in  solutions  that  satisfy  this  constraint  on  the  entire  interval.  The 
optimization  problem  (4.9)  was  solved  with  N  —  50  using  DIDO  version  2003a  [19],  a 
MATLAB  [37]  based  tool  which  implements  the  Legendre  Pseudospectral  collocation 
method,  and  uses  SNOPT  [24]  to  solve  the  resulting  nonlinear  program.  MATLAB 
code  for  solving  this  problem  can  be  found  in  Section  B.l.  Note  that  scaling  factors 
were  used  to  scale  each  variable  in  the  optimization  problem  such  that  the  constraint 
evaluations  and  variable  magnitudes  are  of  similar  orders  of  magnitude.  Values  for 
the  scaling  factors  used  to  solve  this  problem  can  be  found  in  Section  B.l.  The  initial 
guess  for  the  NLP  solver  was  obtained  by  using  a  constant  control  u(t)  =  (0,  0,  0)T 
and  integrating  the  differential  equations  (4.5)  forward  using  ODE45  [37]. 

The  computed  optimal  solutions  are  shown  in  Figures  4. 3-4. 8,  along  with  simula¬ 
tion  results.  The  simulation  results  were  obtained  by  inputting  the  computed  optimal 
control  values  into  a  simulation  which  implements  (4.5)  and  uses  MATLAB’s  ODE45 
for  numerical  integration. 


71 


Momentum  Magnitude  [ft-lbs-sec] 


Momentum  Magnitude  [ft-lbs-sec] 


Figure  4.3:  Simulated,  Optimal  and  Error  Angular  Momentum 
Magnitude  for  Space  Station  Momentum  Dumping  Problem 
with  Continuous  Control,  Using  N=150 

Figure  4.3  shows  the  computed  optimal  objective  function  and  the  simulated  ob¬ 
jective  function  on  the  interval.  It  is  apparent  that  the  computed  results  and  the 
simulated  results  are  in  close  agreement,  the  CMG  momentum  magnitude  was  re¬ 
duced  from  8666  ft-lbs-sec  to  0.1  ft-lbs-sec. 

Figure  4.4  shows  the  computed  optimal  angular  momentum  values  and  the  sim¬ 
ulated  angular  momentum  values  on  the  interval.  Just  as  for  the  magnitude,  the 
angular  momentum  in  each  axis  are  in  close  agreement. 

Figure  4.5  shows  the  optimal  control  computed  by  the  Legendre  Pseudospectral 
method.  As  one  would  expect,  the  control  is  very  nonlinear  which  is  due  to  the 
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Figure  4.4:  Simulated,  Optimal  and  Error  Angular  Momentum  for  Space 
Station  Momentum  Dumping  Problem  with  Continuous 
Control,  Using  N=150 

nonlinear  disturbances  that  are  acting  on  the  spacecraft.  After  the  maneuver  ended, 
at  1800  seconds,  the  control  was  set  to  zero  to  verify  that  a  TEA  was  reached. 

Figure  4.6  depicts  the  optimal  attitude  trajectory  in  Euler  angles.  The  conversion 
from  Euler-Rodrigues  parameters  to  Euler  angles  can  be  found  in  [32].  As  shown  in 
[38]  the  optimal  attitude  trajectory  for  a  desaturation  maneuver  resembles  a  sinusoid 
in  the  roll  axis  and  the  attitude  trajectories  for  the  pitch  and  yaw  axes  are  fairly  flat. 
The  corresponding  angular  rate  trajectories  are  shown  in  Figure  4.7.  As  indicated, 
after  1800  seconds  the  attitude  and  rate  trajectories  are  constant.  This  verifies  that 
a  TEA  was  reached  because  no  control  was  used  to  hold  this  attitude. 
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Figure  4.5:  Computed  Optimal  Control,  Extended  by  Zero  for  t  >  1800, 
for  the  Space  Station  Momentum  Dumping  Problem  with 
Continuous  Control,  Using  N=150. 

Figure  4.8  depicts  the  external  torques  that  are  acting  on  the  spacecraft.  As 
mentioned  earlier,  the  nonlinearity  in  these  torques  account  for  the  nonlinearity  in 
the  optimal  control.  As  further  proof  that  a  TEA  was  reached,  after  the  maneuver 
was  completed  these  torques  are  either  zero  or  counteract  each  other.  This  means 
that  the  right  hand  side  of  (4.1a)  is  zero.  Figure  4.8  shows  that  the  gravity  gradient 
torque  and  the  Euler  torque  combined  were  zero. 

The  optimal  solution  described  in  this  section  is  meaningful  in  the  sense  that 
it  solves  the  problem  of  interest.  However  this  solution  may  not  be  directly  imple- 
mentable  aboard  the  space  station  because  computational  storage  limits  and  pro¬ 
cessing  speed  make  it  difficult  to  command  a  continuous  control  such  as  Figure  4.5. 
Obtaining  implementable  results  requires  that  the  optimal  control  problem  be  solved 
using  piecewise  constant  controls. 
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Figure  4.6:  Simulated,  Optimal  and  Error  Attitude  for  Space  Station 
Momentum  Dumping  Problem  with  Continuous  Control, 
Using  N=150 
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Figure  4.7:  Simulated,  Optimal  and  Error  Angular  Rate  for  Space 
Station  Momentum  Dumping  Problem  with  Continuous 
Control,  Using  N=150 
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Figure  4.8:  Simulated  External  Torques  for  Space  Station  Momentum 
Dumping  Problem  with  Continuous  Control,  Using  N=150 
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4.4  ISS  Momentum  Dumping  Problem  with  Piecewise  Constant 
Control 

In  contrast  to  Section  4.3,  this  version  of  the  space  station  momentum  dumping 
problem  is  posed  such  that  the  control  in  (4.5)  is  a  piecewise  constant  function  on  the 
interval  [£0,£/].  This  is  done  because  it  is  typical  for  the  station  to  perform  a  series 
of  constant  attitude  holds  [39],  which  amount  to  having  a  piecewise  constant  control. 
Again  the  optimal  control  problem  is  be  stated  as 

min  \\h(tf)\\2 

s.t. 


=  Tgg{r) -u(t)x(Ju(t)  +  h(t)) - 

-«(£)>■  £e[£0,£/], 

lr(£) 

—  |(r(£)r(£)T  +  I  +  r(t)x)(u(t) 

-u0(r)),  £e[£0,£/], 

ih(t) 

=  u(t),  t  e  [£0,  £/], 

\m\\2 

—  h"maxi  £  ^  [£o,£/], 

u(t0) 

=  t^O, 

r(£o) 

=  ho, 

h(t0) 

=  h0, 

b(w(tf),r(tf),h(tf )) 

=  0, 

(4.10) 

where  6,  d>o,f0  and  are  given  by  (4.7).  Initial  and  final  conditions  for  the  angular 
rate  and  the  attitude  as  well  as  the  initial  value  for  the  angular  momentum  are 
identical  to  those  in  Section  4.3, 

Since  the  control  is  piecewise  constant  on  the  each  subinterval  the  Legendre 
Pseudospectral  method  with  multiple  subintervals  of  time  must  be  applied  to  dis¬ 
cretize  this  problem.  This  problem  is  posed  on  /  =  5  subintervals  due  to  the  compu¬ 
tational  storage  restrictions  of  onboard  computers.  The  intervals  are  written  as  [to,  £i], 
[£i,£2],  [£2,  £3],  [£3,  £4],  and  [£4,  t5],  with  £0  =  0,£5  =  1800.  Applying  the  pseudospectral 
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direct  transcription  to  (4.8)  results  in  the  NLP 


min  |  [  /i4jv  1 1 2 

s.t. 

2 


N 


i-U-i  J^DjkUik 
k=0 


—  Tgg(rij )  +  hij )  Ui, 


j  =  0,  ...,N,i  =  0,...,4, 

N 

-j—  DjkVik 

3=0 

=  §  irijrij  +  1  +  rij  )  U'u  -  uo  (rij  )  ) 

j  =  0,...,JV,i  =  0,...,4, 

N 

-j—’YhDjkhik 

3=  o 

=  Uii  j  =  0,  =  0,  ...,4, 

(4.11) 


1 1  hi3  1 1 2 

^  hinax  , 

j  —  0, . . . ,  N,i 

Wi-1,N 

—  <Do, 

i  =  !,•••,  4, 

hi-l,N 

hi0 1 

i  =  1,. . .  ,4, 

f i—l,N 

=  Do, 

*  =  !,•••,  4, 

<^0,0 

=  U>0 

O,o 

=  ho 

O 

O 

=  ho 

b(u34tNi  r 4 ,Ni  ^4,Jv) 

=  o, 

where  the  optimization  variables  once  again  are  'jJ13  ,  r,j ,  h%j  and  u%,  j  =  0, . . . ,  N, 
i  =  0,...,4.  The  problem  (4.11)  has  four  new  optimization  variables,  and 

t,4,  because  the  times  for  which  the  control  changes  are  also  to  be  determined.  As  in 
(4.9)  the  path  inequality  constraint  is  only  enforced  at  the  collocation  points  and  it  is 
assumed  that  doing  so  will  result  in  solutions  that  satisfy  this  constraint  on  the  entire 
interval.  The  optimization  problem  (4.11)  was  solved  with  IV  =  30  on  each  subinterval 
using  DIDO  version  2003a  which  relies  on  SNOPT  to  solve  the  resulting  nonlinear 
program.  MATLAB  code  for  solving  this  problem  can  be  found  in  Section  B.2.  Note 
that  scaling  factors  were  used  to  scale  each  variable  in  the  optimization  problem 
such  that  the  constraint  evaluations  and  variable  magnitudes  are  of  similar  orders  of 
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magnitude.  Values  for  the  scaling  factors  used  to  solve  this  problem  can  be  found  in 
Section  B.2. 
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Figure  4.9:  Simulated,  Optimal  and  Error  Angular  Momentum  for 
Space  Station  Momentum  Dumping  Problem  with  Piecewise 
Constant  Control,  Using  N=30  on  Five  Subintervals 


Figure  4.9  shows  the  computed  optimal  objective  function  and  the  simulated  ob¬ 
jective  function  on  the  interval.  As  indicated  the  computed  results  and  the  simulated 
results  are  in  close  agreement,  despite  the  discontinuous  control.  The  momentum 
magnitude  was  reduced  by  almost  6000  ft-lbs-sec,  from  8666  ft-lbs-sec  to  556  ft-lbs- 
sec.  For  this  problem  the  objective  function  value  is  larger  than  the  one  computed  in 
Section  4.3.  Of  course  this  is  to  be  expected  because  the  control  is  restricted  to  be 
from  a  smaller  space  of  functions. 
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Figure  4.10:  Simulated,  Optimal  and  Error  Angular  Momentum  for 
Space  Station  Momentum  Dumping  Problem  with  Piecewise 
Constant  Control,  Using  N=30  on  Five  Subintervals 

Figure  4.10  shows  the  computed  optimal  angular  momentum  values  and  the  sim¬ 
ulated  angular  momentum  values  on  the  interval.  Just  as  for  the  magnitude,  the 
angular  momentum  in  each  axis  are  in  close  agreement.  Figure  4.11  shows  the  opti¬ 
mal  control  computed  by  the  Legendre  Pseudospectral  method. 

Figure  4.12  depicts  the  optimal  attitude  trajectory  in  Euler  angles.  The  attitude 
trajectory  is  similar  to  the  one  in  Section  4.3  but  for  kinks  at  the  control  transitions 
and  the  terminal  values.  As  with  previous  results,  the  roll  trajectory  resembles  a 
sinusoid  while  the  pitch  and  yaw  trajectories  are  fairly  flat.  The  corresponding  an¬ 
gular  rate  trajectories  are  shown  in  Figure  4.13.  Due  to  the  kinks  in  the  attitude  at 
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Figure  4.11:  Computed  Piecewise  Constant  Optimal  Control,  Extended 
by  Zero  for  t  >  1800,  for  Space  Station  Momentum  Dumping 
Problem,  Using  N=30  on  Five  Subintervals. 


the  control  transition  points,  the  rate  is  nearly  discontinuous  at  the  points.  Again, 
after  1800  seconds,  the  control  was  set  to  zero  and  the  attitude  and  rate  trajectories 
remained  constant.  This  verifies  that  a  TEA  was  reached  because  no  control  was  used 
to  hold  this  attitude. 

Figure  4.14  depicts  the  external  torques  that  are  acting  on  the  spacecraft.  As 
further  proof  that  a  TEA  was  reached,  after  the  maneuver  was  completed  these 
torques  are  either  zero  or  counteract  each  other.  Figure  4.14  shows  that  the  combined 
external  torque  goes  to  zero. 
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Figure  4.12:  Simulated,  Optimal  and  Error  Attitude  for  Space  Station 
Momentum  Dumping  Problem  with  Piecewise  Constant 
Control,  Using  N=30  on  Five  Subintervals 
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Figure  4.13:  Simulated,  Optimal  and  Error  Angular  Rate  for  Space 
Station  Momentum  Dumping  Problem  with  Piecewise 
Constant  Control,  Using  N=30  on  Five  Subintervals 
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Figure  4.14:  Simulated  External  Torques  for  Space  Station  Momentum 
Dumping  Problem  with  Piecewise  Constant  Control,  Using 
N=30  on  Five  Subintervals 
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4.4.1  Constraint  on  the  Control  Magnitude 

Some  control  systems  have  a  limit  on  the  amount  of  torque  that  the  attitude  controller 
can  generate  at  any  given  time.  These  limits  manifest  themselves  in  the  form  of  a 
path  inequality  constraint  on  the  control.  Since  the  problem  considered  in  Section  4.4 
uses  piecewise  constant  controls  on  each  subinterval,  constraining  the  control  can  be 
accomplished  by  simple  bounds  on  the  control  on  each  subinterval.  The  constraint 

1 1  fJ'i  1 1 2  ^  ^ max  >  i  1,...,5  (4.12) 

was  added  to  the  NLP  (4.11)  to  produce  the  results  shown  in  Figures  4.15  through 
4.20  for  umax  =  200  ft-lbs.  MATLAB  code  for  solving  this  problem  can  be  found  in 
Section  B.2.  Again,  note  that  scaling  factors  were  used  to  scale  each  variable  in  the 
optimization  problem  such  that  the  constraint  evaluations  and  variable  magnitudes 
are  of  similar  orders  of  magnitude.  Values  for  the  scaling  factors  used  to  solve  this 
problem  can  be  found  in  Section  B.l. 
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Figure  4.15:  Simulated,  Optimal  and  Error  Angular  Momentum 
Magnitude  for  Space  Station  Momentum  Dumping  Problem 
with  Constrained  Piecewise  Constant  Control,  Using  N=30 
on  Five  Subintervals 


Simulated 
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Figure  4.16:  Simulated,  Optimal  and  Error  Angular  Momentum 
for  Space  Station  Momentum  Dumping  Problem  with 
Constrained  Piecewise  Constant  Control,  Using  N=30  on 
Five  Subintervals 
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Figure  4.17:  Computed  Piecewise  Constant  Optimal  Control,  Extended 
by  Zero  for  t  >  1800,  for  Space  Station  Momentum  Dumping 
Problem,  Using  N=30  on  Five  Subintervals 
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Figure  4.18:  Simulated,  Optimal  and  Error  Attitude  for  Space  Station 
Momentum  Dumping  Problem  with  Constrained  Piecewise 
Constant  Control,  Using  N=30  on  Five  Subintervals 


Simulated 


90 


Rate  [deg/sec]  Rate  [deg/sec] 


Time  [sec] 


Figure  4.19:  Simulated,  Optimal  and  Error  Angular  Rate  for  Space 
Station  Momentum  Dumping  Problem  with  Constrained 
Piecewise  Constant  Control,  Using  N=30  on  Five 
Subintervals 
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Figure  4.20:  Simulated  External  Torques  for  Space  Station  Momentum 
Dumping  Problem  with  Constrained  Piecewise  Constant 
Control,  Using  N=30  on  Five  Subintervals 
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4.5  ISS  Momentum  Dumping  Problem  with  Control  Law 

This  problem  is  a  variation  of  the  one  described  in  Section  4.3  but  with  a  control  law 
introduced  into  the  dynamics.  Using  a  control  law  rather  than  a  torque  command, 
as  done  in  Sections  4. 3-4. 4,  allows  the  ISS  to  be  controlled  by  an  attitude  command. 
Typically  this  control  law  is  used  to  drive  the  spacecraft  to  some  desired  attitude  by 
inputting  an  attitude  command  into  the  controller.  This  is  an  important  variation 
to  the  momentum  dumping  problem  because  the  ISS  actually  controls  its  attitude 
through  this  type  of  control  law.  Therefore,  obtaining  an  optimal  attitude  command, 
as  opposed  to  an  optimal  torque  control,  has  more  practical  value.  The  control 
law  used  here  is  a  proportional  derivative  control  law  which  is  described  in  [51]. 
This  control  law  is  formulated  in  terms  of  quaternions.  The  relationship  between 
quaternions  and  Euler-Rodriguez  parameters,  which  are  used  here,  can  be  found  in 
[32], 

The  state  variables  are  the  angular  rate  u,  the  attitude  r  and  the  CMG  angular 
momentum  h.  The  control  variable  here  is  rcj  instead  of  u  due  to  the  fact  that  the 
control  for  this  problem  is  actually  the  desired  attitude  not  the  time  derivative  of 
momentum  as  in  Section  4.3. 

min  \\h{tf)\\2 
s.t. 


JfMt) 

=  Tgg(r)  ~uj(t)x(J<jj(t)  +  h(t))  -  fth(t ),  t  G  [t0,tf], 

iAt) 

=  l(r(t)r(t)T  +  I  +  r(t)x)(u(t) -w0(r)),  te[t0,tf}, 

ftA) 

=  k\Jioerr  (cu,  r d)  +  k2Jrerr(r,rd),  t  G  [t0,tf], 

\m\\* 

5s  hmaxi  f  £  [to,tf], 

u(t0) 

=  ^0, 

r(to) 

=  h0, 

h(t0) 

—  ^0; 

b(ui(tf),r(tf),h(tf))  =  0, 


(4.13) 
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where  :  M  i— >  M3  is  the  attitude  command  and 
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c3 

c2 

T 

R 


=  /  + 


2  /.j.Xy.x  _  rx\ 


1  +rTr 


C  (o,0,l)T, 
C  (o,i,o)T, 


=  C*(JC3), 


Ro  = 


R  (0,0,l)T, 


err  ^  H-  ^orbR/2i 


D  =  CRT , 


f  err  ~  |(1  +  -Dll  +  D22  +  D33)  l^2 


(  \ 

d23  —  d32 

D31  —  d13 


D\2  —  D  21 


and 


LU0  ^orbC'2- 


(4.14) 


In  (4.14),  R2  is  the  second  column  of  the  rotation  matrix  associated  with  the  attitude 
control  r<2 .  Again,  b,  ou0,  fo  and  h o  are  given  by  (4.7).  The  gains  for  the  CMG  controller 
k\  and  k2  are  given  as 

h  =  0.0632, 
k2  =  0.002. 

The  control  law  described  here  forces  the  spacecraft  to  go  to  the  desired  attitude 
with  a  fixed  angular  rate.  This  is  a  restrictive  model  because  real  systems  tend  to 
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reach  the  required  rate  only  at  the  end  of  the  maneuver.  However,  this  control  law  is 
consistent  with  the  one  used  in  [38]  to  control  the  same  system. 

MATLAB  code  for  solving  this  problem  can  be  found  in  Section  B.3.  Note  that 
scaling  factors  were  used  to  scale  each  variable  in  the  optimization  problem  such  that 
the  constraint  evaluations  and  variable  magnitudes  are  of  similar  orders  of  magnitude. 
Values  for  the  scaling  factors  used  to  solve  this  problem  can  be  found  in  Section  B.3. 
Results  for  this  problem  for  N  =  150  are  shown  in  Figures  4.21  through  4.27. 
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Figure  4.21:  Simulated,  Optimal  and  Error  Angular  Momentum 
Magnitude  for  Space  Station  Momentum  Dumping  Problem 
with  Control  Law,  Using  N=150 
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Figure  4.22:  Simulated,  Optimal  and  Error  Angular  Momentum  for  Space 
Station  Momentum  Dumping  Problem  with  Control  Law, 
Using  N=150 
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Figure  4.23:  Simulated,  Optimal  and  Error  Angular  Momentum  Time 
Derivative  for  Space  Station  Momentum  Dumping  Problem 
with  Control  Law,  Using  N=150 
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Figure  4.24:  Computed  Optimal  Attitude  Command  Control,  Extended 
by  Zero  for  t  >  1800,  for  the  Space  Station  Momentum 
Dumping  Problem,  Using  N=150 
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Figure  4.25:  Simulated,  Optimal  and  Error  Attitude  for  Space  Station 
Momentum  Dumping  Problem  with  Control  Law,  Using 


N=150 
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Figure  4.26:  Simulated,  Optimal  and  Error  Angular  Rate  for  Space 
Station  Momentum  Dumping  Problem  with  Control  Law, 
Using  N=150 
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Figure  4.27:  Simulated  External  Torques  for  Space  Station  Momentum 
Dumping  Problem  with  Control  Law,  Using  N=150 
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The  results  presented  in  this  chapter  represent  a  solution  which  is  accurate  in  the 
sense  that  it  has  been  verified  through  simulation,  and  practical  in  the  sense  that 
piecewise  constant  controls  are  implementable  aboard  the  space  station.  This  chap¬ 
ter  demonstrated  the  utility  of  the  Legendre  Psendospectral  method  for  solving  an 
optimal  control  problem.  Dne  to  the  accuracy  properties  of  this  method,  this  problem 
was  solved  using  few  collocation  points.  The  pseudospectral  collocation  method  was 
also  versatile  enough  to  account  for  time  dependent  and  piecewise  constant  control 
problems  with  little  modification  to  the  direct  transcription  of  the  optimal  control 
problem. 
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Chapter  5 

Conclusions  and  Future  Work 

In  this  thesis  pseudospectral  collocation  methods  for  the  direct  transcription  of  op¬ 
timal  control  problems  were  presented.  It  was  shown  that  these  methods  exhibit 
properties  which  make  it  possible  to  relate  the  discretized  NLP  to  the  infinite  di¬ 
mensional  OCP.  This  was  done  by  constructing  a  new  adjoint  mapping  that  relates 
the  Lagrange  multipliers  in  the  discretized  optimal  control  problem  to  the  adjoint 
variables  corresponding  to  the  infinite  dimensional  optimal  control  problem. 

Through  this  adjoint  estimation  procedure,  error  estimates  between  the  computed 
solution  and  the  true  solution  to  the  optimal  control  problem  were  derived  for  linear- 
quadratic  optimal  control  problems.  It  was  shown  that  the  conditioning  of  the  opti¬ 
mality  system  matrix  plays  an  important  role  in  obtaining  accurate  solutions. 

These  methods  were  applied  to  the  International  Space  Station  momentum  dump¬ 
ing  problem  to  demonstrate  their  utility  for  difficult  problems.  The  solutions  obtained 
here  are  of  great  practical  value  as  they  may  be  directly  implementable  aboard  oper¬ 
ational  spacecraft. 

One  suggestion  for  future  work  is  to  compute  error  estimates  for  the  optimal  con¬ 
trol  for  nonlinear  optimal  control  problems.  This  would  be  an  significant  result  in  the 
sense  that  many  important  applications  are  using  pseudospectral  collocation  meth¬ 
ods  to  solve  optimal  control  problems.  Such  error  estimate  would  give  the  scientists, 
engineers  and  mathematicians  who  use  these  methods  more  confidence  in  computed 
solutions. 

Another  suggestion  for  future  work  would  be  to  develop  a  more  robust  optimal 
control  solver  that  implements  pseudospectral  methods.  Current  packages  either  do 
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not  give  the  user  sufficient  insight  into  the  nonlinear  programming  solver,  do  not 
allow  the  user  to  supply  analytical  derivatives  or  do  not  allow  the  user  to  choose 
which  pseudospectral  collocation  method  is  used. 

A  final  suggestion  would  be  to  use  pseudospectral  methods,  or  another  direct  tran¬ 
scription  method  for  that  matter,  to  solve  extensions  of  the  ISS  momentum  dumping 
problem,  as  solutions  to  this  problem  can  be  quite  useful  in  the  aerospace  industry. 
Extensions  to  the  ISS  momentum  dumping  problem  include  solving  for  a  discrete 
(piece-wise  constant)  attitude  command  control  that  will  remove  built-up  momen¬ 
tum,  solving  for  an  optimal  control  that  is  robust  to  unmodeled  dynamics  or  changes 
in  the  boundary  conditions,  and  attitude  maneuvers  during  payload  operations. 
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Appendix  A 


DIDO:  A  Tool  for  Direct  and  InDirect 
Optimization 


DIDO  [19]  is  a  MATLAB  [37]  based  tool  which  implements  the  Legendre  Pseudospectral 
method  described  in  Section  3.  The  resulting  nonlinear  programming  problem  is 
solved  using  SNOPT  [24]  which  is  interfaced  through  MATLAB  via  an  optimization 
tool  called  TOMLAB  [31].  TOMLAB  implements  various  optimization  solvers,  one 
of  which  is  SNOPT,  using  MEX-hles  which  call  supporting  Fortran  routines.  DIDO 
discretizes  the  infinite  dimensional  optimal  control  problem  and  sends  the  resulting 
NLP  into  TOMLAB  to  obtain  solutions. 

DIDO  is  set  up  to  solve  optimal  control  problems  of  the  following  form 

mmE(y(t0),y(tf),t0,tf)+  f  F(y(t),u(t),t)dt ,  (A.la) 

Jto 

subject  to  the  dynamic  constraints 

=  f(y(t),u(t),t),  (A. lb) 

the  event  constraints 

ei  <  e(y(tQ),y(tf),t0,  tf)  <  eu,  (A.lc) 

the  path  inequality  constraints 

hi  <  h(y(t),u(t),t)  <  hu,  (A. Id) 

and  the  state  and  control  bounds 

yi  <  y(t )  <  yu, 

Ui  <  u(t )  <  Uu. 


(A.le) 


Ill 


In  (A.l)  we  have  used  the  notation  applied  in  [19],  which  is  slightly  different  from 
the  one  used  in  previous  chapters.  To  use  this  tool,  one  must  define  at  least  two,  up 
to  four,  auxiliary  functions.  An  M-file  which  defines  the  objective  function  (A. la). 
This  file  must,  given  state,  control  and  time  values,  return  the  end-time  cost  function 
value  E  and  the  integral  cost  function  value  F  at  each  collocation  point.  In  addition, 
an  M-file  which  defines  the  dynamics  function  must  be  defined.  This  file  must,  given 
state,  control,  time  and  state  time  derivative  values,  return  the  difference  between  the 
approximated  state  time  derivative  yN  and  the  right  hand  side  function  /  at  each 
collocation  point.  The  optional  event  function  can  be  defined,  which,  given  state  and 
time  initial  and  final  values,  return  the  value  of  the  event  function  e.  The  optional 
path  inequality  function  can  be  defined,  which,  given  state,  control  and  time  values, 
return  the  value  of  the  path  function  h.  Bounds  on  each  function,  the  states  and  the 
controls  are  passed  directly  to  the  DIDO  function  call.  Examples  of  each  function 
can  be  found  in  Section  B. 

DIDO  is  a  relatively  easy  tool  to  use  because  the  setup  time  involved  in  solving 
a  given  optimal  control  problem  can  be  small  as  compared  to  most  alternatives.  The 
coding  that  is  required  of  a  user  is  minimal  in  the  sense  that  only  a  few,  possibly 
small,  M-files  need  to  be  programmed.  Additionally,  the  user  interface,  or  DIDO 
function  call,  is  very  straight  forward  and  resembles  that  of  any  standard  MATLAB 
function  call. 

Despite  these  nice  characteristics,  there  are  some  unattractive  aspects  to  DIDO. 
The  first  is  that,  as  the  reader  may  have  noticed,  none  of  the  user  supplied  functions 
output  derivative  information.  Not  only  do  they  not  require  derivative  information, 
the  tool  is  not  set  up  to  allow  the  use  of  user  supplied  derivatives.  Instead,  finite 
difference  derivatives  are  computed  within  SNOPT.  This  has  several  consequences. 
First  and  foremost,  it  is  well  known  that  the  performance  of  an  optimization  algo¬ 
rithm  can  be  severely  slowed  down  when  finite  difference  derivatives  are  used.  As  a 
consequence,  solving  relatively  small  problems  can  be  very  time  consuming  ordeal, 
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even  when  exact  derivatives  may  be  easy  to  compute.  Secondly,  the  accuracy  of 
solutions  are  impacted  because  the  quality  of  derivative  approximations  are  limited 
by  the  error  in  function  values.  In  our  context,  error  in  function  values  are  related 
to  machine  precision  e.  For  example,  the  error  in  derivative  approximations,  when 
forward  finite  differences  are  used,  is  on  the  order  of  ^/e  [28].  Therefore,  solutions 
that  are  obtained  by  using  DIDO  may  be  less  accurate  than  expected,  especially  if 
high  order  pseudospectral  collocation  discretizations  are  used. 

Another  disadvantage  is  that  DIDO  does  not  give  the  user  the  opportunity  to 
control  the  NLP  solver.  Stopping  tolerances,  iteration  limits  and  the  like  are  all 
set  within  DIDO.  As  a  consequence,  accuracy  expected  from  high  order  Legendre 
Pseudospectral  collocation  methods  may  be  polluted  by  coarse  optimization  stopping 
tolerances.  The  user,  while  free  to  select  the  degree  N  of  the  Legendre  Pseudospectral 
has  no  opportunity  to  adjust  the  accuracy  of  the  NLP  solve  to  match  the  expected 
discretization  error. 
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Appendix  B 


MATLAB  Code  for  the  Solution  of  the  ISS 
Momentum  Dumping  Problem  Using  DIDO 


This  Appendix  contains  the  MATLAB  code  that  may  be  used  to  produce  the  results 
for  the  ISS  momentum  dumping  problem  as  solved  by  DIDO  in  Section  4. 


B.l  ISS  Momentum  Dumping  Problem  with  Continuous  Control 

The  main  hie  which  calls  DIDO  is  shown  below. 


% 

°/0  issMain 

l 

°/0  solve  the  ISS  momentum  dumping  problem  for 
°/0  continuous  control 

l 

1  min  | I h(tf ) |  | 

1 

1  s.t. 

7.  Jw’  =  Tgg(r)  -  w  x  (Jw+h)  -  u 
%  h’  =  u 

°/o  r’  =  0.5  *  (  rr’  +  I  +  skew(r))(w  -  w(r)) 
%  I  I h | |  <=  10000 
% 


% - 

°/„  Set  global  variables 

% - 

global  ws  hs  rs  us  w_orb  J  iJ 

°/„  Scaling  factors 
ws  =  le-3; 
hs  =  le3; 
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rs  =  leO; 
us  =  leO; 


°/„  Orbital  rate 

w_orb  =  0 . 06511* (pi/180) ; 


7.  Inertia  matrix 

J  =  [2 . 807019116160000e+007  4 . 822509936000001e+005  -1 . 716750944480000e+007 
4 . 82250993600000 le+005  9 . 514463934400001e+007  6 . 026044480000001e+004 

-1 . 7 16750944480000e+007  6 . 026044480000001e+004  7 . 659440133600001e+007] ; 

7„  Inertia  matrix  inverse 

iJ  =  [4. 128859554604031e-008  -2 . 151370870617538e-010  9 . 254401181964391e-009 

-2 . 15 13708706 17538e-0 10  1 . 051143985685158e-008  -5 . 648966425550106e-011 

9 . 25440 118 196439 le-009  -5 . 648966425550106e-011  1 . 513006699674998e-008] 


7. - 

7„  Set  input  functions 

% - 

iss.cost  =  ’ issCost’; 
iss. dynamics  =  ’ issDynamics ’ ; 
iss.path  =  ’ issPath’; 
iss. events  =  ’ issEvents ’ ; 


°/o - 

°/„  Set  variable  parameters 

y. - 

to  =  0; 

tf  =  1800; 

wO  =  [-9 . 5380685844896e-006  -1 . 1363312657036e-003  5 . 3472801108427e-006] /ws ; 
wlb  =  inf*[-l  -1  -l]/ws; 
wub  =  inf*[l  1  l]/ws; 
wf  =  wO ; 

hO  =  [5e3  5e3  5e3]/hs; 
hlb  =  [-le4  -le4  -le4]/hs; 
hub  =  [le4  le4  le4]/hs; 

rO  =  [2.9963689649816e-003  1 . 5334477761054e-001  3 . 8359805613992e-003] /rs ; 
rib  =  inf*[-lelO  -lelO  -lelO]/rs; 
rub  =  inf*[lelO  lelO  lelO]/rs; 
rf  =  rO; 

% - 

7.  Set  time  bounds 

% - 

knots . locations  =  [tO  tf]  ; 
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knots . definitions  =  {’hard’,  ’hard’}; 
knots . bounds . lower  =  [tO  tf] ; 
knots . bounds .upper  =  [tO  tf]  ; 
knots . numNodes  =  [150]; 


7. - 

7.  Set  variable  bounds 

% - 

bounds . lower . states ( : , 1) 
bounds .upper . states ( : , 1) 
bounds . lower . states ( : , 2) 
bounds . upper . states ( : , 2) 
bounds . lower . states ( : , 3) 
bounds . upper . states ( : , 3) 


[wO  hO  rO] ’ ; 

[wO  hO  rO] ’ ; 
[wlb  hlb  rib] ’ ; 
[wub  hub  rub] ’ ; 
[wlb  hlb  rib] ’ ; 
[wub  hub  rub] ’ ; 


bounds . lower . controls  =  -inf  *  [111]’; 
bounds .upper . controls  =  inf  *  [111]’; 


% - 

7.  Set  constraint  bounds 

7. - 

bounds . lower .path  =  [0] ’ ; 
bounds .upper .path  =  [le8  /  hs~2] ’ ; 
bounds . lower . events  =0*  [111111]’; 
bounds .upper . events  =0*  [111111]’; 


7. - 

7.  Provide  a  guess 

7. - 

load  iss_cont_guess ; 


7. - 

7.  Call  DIDO 

7. - 

[cost,  primal, dual]  =  dido(iss,  knots,  bounds,  guess); 


The  main  file  sets  the  dynamics  function,  the  events  function,  and  the  path  func¬ 
tion  for  the  DIDO  tool.  It  also  defines  global  variables,  sets  bounds  on  variables 
and  function  right-hand  sides,  and  loads  an  initial  guess  from  a  data  file.  Note  that 
weighting  parameters  were  used  to  rescale  each  variable  to  be  on  the  same  order  of 
magnitude.  This  procedure  is  discussed  in  [7,  19]. 
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The  function  which  implements  the  differential  equation  is  shown  below. 


function  [R]  =  issDynamics (primal) 

7. 

°/„  ISS  momentum  control  dynamics 

l 

°/„  Jw’  =  Tgg(r)  -  w  x  (Jw+h)  -  u 
%  h’  =  u 

°/„  r’  =  0.5  *  (  rr’  +  I  +  skew(r))(w  -  w(r)) 

°/„  global  variables 

global  ws  hs  rs  us  w_orb  J  iJ 

°/„  initialize  variables 
w  =  primal . states (1 : 3 ,: )  *  ws; 
h  =  primal . states (4 : 6 ,: )  *  hs; 
r  =  primal . states (7 : 9 ,: )  *  rs; 
u  =  primal . controls  *  us; 

N  =  length (w( 1 , : ) ) ; 

R  =  zeros(9,N) ; 

for  i  =  1:N 

l  compute  auxiliary  values 

C  =  eye (3)  +  2/(l+r ( : , i) 5 *r ( : , i) ) * (skew(r ( : , i) ) * . . . 

skew(r(:,i))  -  skew(r ( : , i) ) ) ; 

C2  =  C( : ,2) ; 

C3  =  C( : ,3) ; 

wp  =  -w_orb  *  C2; 

Tgg  =  3*w_orb~2* (cross (C3 , ( J*C3) ) ) ; 

%  compute  differential  constraint 

R(:,i)  =  primal . statedots (:, i)  -  [  iJ*(Tgg  -  ... 

cross (w(: ,  i) ,J*w(: , i)+h( : ,  i) )  -  u(:,i))/ws; 

(u( : ,i))/hs; 

(0 . 5*(r ( : , i) *r ( : , i) ’  +  eye(3)  +  skew(r ( : , i) ) )*(w( : , i)-wp) )/rs] ; 

end 


As  discussed  in  [19]  it  is  more  computationally  efficient  to  vectorize  all  functions. 
In  order  to  do  this,  the  above  dynamics  function  must  be  expanded  out  component¬ 
wise.  For  the  sake  of  brevity,  this  exercise  is  left  out  of  this  thesis. 
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The  function  which  enforces  the  final  time  constraint  is  shown  below. 


function  [E]  =  issEvents (primal) 

% 

°/0  this  function  enforces  that  a  TEA  is  reached 

7.  at  the  final  time  by  requireing  that  the  ODE 

°/„  is  zero  at  the  final  time  when  the  control  is  zero 

% 


7.  global  variables 

global  ws  hs  rs  us  w_orb  J  iJ 

7.  initialize  variables 
w  =  primal . states (1 : 3 , end)  *  ws; 
h  =  primal . states (4 : 6 , end)  *  hs; 
r  =  primal . states (7 : 9 , end)  *  rs; 

7.  compute  auxiliary  values 


c 

=  eye (3)  + 

2/ (l+r’*r)*(skew(r)*skew(r)  -  skew(r)); 

C2 

=  C(: ,2) ; 

C3 

=  C(: ,3); 

wp 

=  -w_orb  * 

C2 ; 

Tgg  =  3*w_orb~2*(cross(C3, (J*C3))) ; 

7.  compute  differential  constraint 
E  =  [  iJ*(Tgg  -  cross(w, J*w+h))/ws; 

(0.5*(r*r’  +  eye(3)  +  skew(r))*(w-wp))/rs] ; 


The  above  function  sets  the  constraint  that  a  toque  equilibrium  attitude  must  be 
reached  at  the  final  time.  As  described  in  Section  4.2,  this  is  equivalent  to  forcing 
the  right-hand  side  of  the  differential  equation  to  zero. 

The  function  that  implements  the  path  inequality  constraint  is  shown  below. 


function  [P]  =  issPath(primal) 

7. 

7.  this  function  enforces  the  path  inequality  constraint 
7.  that  |  |  h  |  I  is  less  than  or  equal  to  10000  at  each  time 

7. 
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°/„  global  variables 

global  ws  hs  rs  us  w_orb  J  iJ 

°/„  2-norm-squared  of  h  at  each  time 

P  =  [(primal . states (4, : )*hs) . ~2  +  (primal . states (5 ,:) *hs) . "2  +. . . 
(primal . states (6 ,:) *hs) . ~2]  /  hs~2; 


The  above  function  enforces  the  path  inequality  constraint  on  the  norm  of  the 
angular  momentum. 

The  function  that  implements  the  final  time  cost  function  on  the  norm  of  the 
angular  momentum  is  shown  below. 


function  [M,  I]  =  issCost (primal) 

7. 

7.  this  function  is  the  final  time  objective  function 
7.  on  the  2-norm  of  h  (squared) 

7. 


7.  global  variables 

global  ws  hs  rs  us  w_orb  J  iJ 

7.  final  time  function 

M  =  hs~2  *  primal . states (4 : 6 , end) ’ *primal . states (4 : 6 , end) ; 

7.  no  integral  function 
I  =  zeros(size(  primal . states (1 ,: )  )  ); 


B.2  ISS  Momentum  Dumping  Problem  with  Piecewise  Constant 
Control 

The  main  file  which  calls  DIDO  is  shown  below. 


7. 

7.  issMain 

7. 
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°/„  solve  the  ISS  momentum  dumping  problem  for 
°/„  pw  constant  control 

l 

1  min  | I h(tf ) |  | 

1 

1  s.t. 

7.  Jw’  =  Tgg(r)  -  w  x  (Jw+h)  -  u 
%  h’  =  u 

°/„  r’  =  0.5  *  (  rrJ  +  I  +  skew(r))(w  -  w(r)) 
l  | | h | |  <=  10000 

°/o 


°/o - 

°/„  Set  control  bound 
°/„  0  =  no 
%  1  =  yes 

% - 

bounded_control  =  0; 


°/o - 

°/„  Set  global  variables 
% - 

global  ws  hs  rs  us  N1  N2  N3  N4  N5  J  iJ  w_orb 

°/„  N  on  each  subinterval 
Ml  =30; 

M2  =30; 

M3  =30; 

M4  =  30 ; 

M5  =30; 

°/„  Scaling  factors 
ws  =  le-2; 
hs  =  le3; 
rs  =  le-2; 
us  =  leO; 

°/„  Orbital  rate 

w_orb  =  0 . 06511* (pi/180) ; 

°/„  Inertia  matrix 

J  =  [2 . 807019116160000e+007  4 . 822509936000001e+005  - 
4 . 82250993600000 le+005  9 . 514463934400001e+007 
-1 . 716750944480000e+007  6 . 026044480000001e+004 


1 . 716750944480000e+007 
6 . 02604448000000 le+004 
7 . 659440133600001e+007] ; 
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7.  Inertia  matrix  inverse 


iJ  = 


[4. 128859554604031e-008 
-2. 151370870617538e-010 
9.254401181964391e-009 


-2 . 15 13708706 17538e-0 10 
1 . 051 143985685 158e-008 
-5 . 648966425550 106e-0 11 


9 . 25440 118 196439 le-009 
-5 . 648966425550 106e-0 11 
1.513006699674998e-008] ; 


7. - 

7„  Set  input  functions 

7. - 

iss.cost  =  ’issCost’; 
iss. dynamics  =  ’ issDynamics ’ ; 
iss.path  =  ’issPath’; 
iss. events  =  ’ issEvents ’ ; 


7. - 

7„  Set  variable  parameters 

7. - 

tO  =  0; 
tf  =  1800; 

wO  =  [-9 . 5380685844896e-006  -1 . 1363312657036e-003  5 . 3472801108427e-006] /ws ; 
wlb  =  [-1  -1  -l]/ws; 
wub  =  [1  1  l]/ws; 
wf  =  wO ; 

hO  =  [5e3  5e3  5e3]/hs; 
hlb  =  [-le4  -le4  -le4]/hs; 
hub  =  [le4  le4  le4]/hs; 

rO  =  [2.9963689649816e-003  1 . 5334477761054e-001  3 . 8359805613992e-003] /rs ; 
rib  =  [-lelO  -lelO  -lelO]/rs; 
rub  =  [lelO  lelO  lelO]/rs; 
rf  =  rO; 

% - 

7.  Set  time  bounds 

% - 

small  =  100*eps; 

knots . locations  =  [tO  360  720  1080  1440  tf]  ; 

knots . definitions  =  {’hard’,  ’soft’,  ’soft’,  ’soft’ , ’soft’ , ’hard’}; 
knots . bounds . lower  =  [tO  tO+small  540  1000+small  1400  tf] ; 
knots . bounds .upper  =  [tO  540-small  1000  1400-small  tf-small  tf] ; 
knots . numNodes  =  [N 1  N2  N3  N4  N5] ; 


°/o - 

7.  Set  variable  bounds 

% - 

bounds . lower . states (:, 1)  =  [wO  hO  rO]  ’ ; 
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bounds . upper . states ( : 
bounds . lower . states ( : 
bounds . upper . states ( : 
bounds . lower . states ( : 
bounds . upper . states ( : 


,1)  = 

[wO  hO  rO] ’ ; 

,2)  = 

[wlb 

hlb 

rib] 

,2)  = 

[wub 

hub 

rub] 

,3)  = 

[wlb 

hlb 

rib] 

,3)  = 

[wub 

hub 

rub] 

7.  note:  parameters  act  as  controls  here, 

7„  but  the  must  be  a  control  variable 
7„  "placeholder"  which  is  not  iterated  upon 
if  bounded_control  ==  0 

bounds . lower . parameters 
bounds . upper . parameters 
elseif  bounded_control  ==  1 
bounds . lower . parameters 
bounds . upper . parameters 

end 


-lelO  *  ones(5*3,l); 
lelO  *  ones(5*3,l); 

-200  *  ones(5*3,l); 
200  *  ones(5*3, 1) ; 


bounds . lower . controls  =  0; 
bounds .upper . controls  =  0; 


°/o - 

°/„  Set  constraint  bounds 

% - 

bounds . lower .path  =  [0]’; 

bounds .upper .path  =  [le2]  ’ ; 

bounds . lower . events  =0*  [111111]’; 

bounds .upper . events  =0*  [111111]’; 


°/o - 

°/„  Provide  a  guess 

°/o - 

load  iss_dis_guess ; 


% - 

%  Call  DIDO 

% - 

[cost,  primal, dual]  =  dido(iss,  knots,  bounds,  guess); 


This  main  file  does  the  same  things  as  the  one  shown  in  Section  B.l  except  for 
two  main  things.  The  first  is  that  the  time  interval  is  broken  up  into  five  subintervals 
so  that  a  constant  control  can  be  used  on  each  subinterval.  This  leads  to  the  second 
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difference.  In  terms  of  the  DIDO  tool,  the  control  in  this  setting  is  not  used  and 
the  variable  acting  as  the  control  is  viewed  as  a  parameter.  DIDO  must  evaluate  the 
control  at  every  collocation  point  and  it  forces  these  variable  to  be  variable  in  the 
optimization,  eventhough  the  user  may  want  the  value  to  remain  constant  over  each 
subinterval.  The  use  of  parameters  alleviates  this  problem.  Accomplishing  the  bound 
on  the  control  can  be  accomplished  by  changing  the  parameter  bound  in  the  main 
driver  hie. 

The  function  which  implements  the  differential  equation  is  shown  below. 


function  [R]  =  issDynamics (primal) 

7. 

7.  ISS  momentum  control  dynamics 

l 

°/„  Jw’  =  Tgg(r)  -  w  x  (Jw+h)  -  u 
%  h’  =  u 

°/0  r’  =  0.5  *  (  rr’  +  I  +  skew(r))(w  -  w(r)) 

7.  global  variables 

global  ws  hs  rs  us  N1  N2  N3  N4  N5  J  iJ  w_orb 

7.  initialize  variables 

w  =  primal . states (1 : 3 ,: )  *  ws; 

h  =  primal . states (4 : 6 )  *  hs; 

r  =  primal . states (7 : 9 ,: )  *  rs; 

u  =  reshape (primal .parameters , 3, 5)  *  us; 

u  =  [kron(u( : , 1) , ones (1 ,N1) )  kron(u( : , 2) , ones (1 ,N2) )  ... 

kron(u( : , 3) , ones (1 ,N3) )  kron(u( : ,4) , ones (1 ,N4) )  ... 
kron(u( : , 5) , ones (1 ,N5) )] ; 

N  =  length (w( 1 , : ) ) ; 

R  =  zeros(9,N) ; 

for  i  =  1:N 

%  compute  auxiliary  values 
C  =  eye(3)  +  2/(l+r ( : , i) J *r ( : , i) )* .  .  . 

(skew(r ( : , i) ) *skew(r ( : , i) )  -  skew(r( : ,i))) ; 

C2  =  C(: ,2) ; 

C3  =  C(: ,3) ; 

wp  =  -w_orb  *  C2; 

Tgg  =  3*w_orb~2*(cross(C3, (J*C3))) ; 
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%  compute  differential  constraint 

R(:,i)  =  primal . statedots (:, i)  -  [  iJ*(Tgg  -  ... 

cross (w( : ,  i) , J*w( : , i)+h( : ,i))  -  u(:,i))/ws; 

(u(: ,i))/hs; 

(0.5*(r(: ,i)*r(: ,i) ’  +  eye(3)  +  skew(r ( : , i) ) ) * (w( : , i)-wp) )/rs] ; 

end 


The  events  function,  the  path  function  and  the  cost  function  for  this  problem  are 
the  same  as  the  ones  stated  in  Section  B.l. 


B.3  ISS  Momentum  Dumping  Problem  with  Control  Law 

The  main  file  which  calls  DIDO  is  shown  below. 

% 

°/„  issMain 

°/o 

7.  solve  the  ISS  momentum  dumping  problem  for 
°/„  control  law 

°/o 

°/o  min  |  I h(tf )  |  | 

l 

1  s.t. 

7.  Jw’  =  Tgg(r)  -  w  x  (Jw+h)  -  u 

7.  h’  =  kl*J*w_e  +  k2*J*r_e 

%  r’  =0.5*(rr’  +1+  skew(r))(w  -  w(r)) 

7.  I  I  h  |  |  <=  10000 

l 


7. - 

7.  Set  global  variables 

l - 

global  ws  hs  rs  us  kl  k2  J  iJ  w_orb 

7.  Scaling  factors 
ws  =  le-2; 
hs  =  le3; 
rs  =  le-1; 
us  =  le-1; 


7.  Orbital  rate 
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w_orb  =  0 . 06511* (pi/180) ; 


°/„  Inertia  matrix 

J  =  [2 . 807019116160000e+007  4 . 822509936000001e+005  -1 . 716750944480000e+007 
4 . 82250993600000 le+005  9 . 514463934400001e+007  6 . 026044480000001e+004 

-1 . 716750944480000e+007  6 . 026044480000001e+004  7 . 659440133600001e+007] ; 

7.  Inertia  matrix  inverse 

iJ  =  [4. 128859554604031e-008  -2 . 151370870617538e-010  9 . 254401181964391e-009 

-2 . 15 13708706 17538e-0 10  1 . 051143985685158e-008  -5 . 648966425550106e-011 

9 . 25440 118 196439 le-009  -5 . 648966425550106e-011  1 . 513006699674998e-008] ; 


7„  Controller  gains 
kl  =  0.0632; 
k2  =0.002; 


7. - 

7„  Set  input  functions 

% - 

iss.cost  =  ’ issCost’; 
iss. dynamics  =  ’ issDynamics ’ ; 
iss.path  =  ’ issPath’; 
iss. events  =  ’ issEvents ’ ; 


°/o - 

°/„  Set  variable  parameters 

y. - 

to  =  0; 

tf  =  1800; 

wO  =  [-9 . 5380685844896e-006  -1 . 1363312657036e-003  5 . 3472801108427e-006] /ws ; 
wlb  =  [-1  -1  -l]/ws; 
wub  =  [1  1  l]/ws; 
wf  =  wO ; 

hO  =  [5e3  5e3  5e3]/hs; 
hlb  =  [-le4  -le4  -le4]/hs; 
hub  =  [le4  le4  le4]/hs; 

rO  =  [2.9963689649816e-003  1 . 5334477761054e-001  3 . 8359805613992e-003] /rs ; 
rib  =  [-lelO  -lelO  -lelO]/rs; 
rub  =  [lelO  lelO  lelO]/rs; 
rf  =  rO; 

°/o - 

y.  Set  time  bounds 

% - 

knots . locations  =  [tO  tf]  ; 
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knots . definitions  =  {’hard’,  ’hard’}; 
knots . bounds . lower  =  [tO  tf] ; 
knots . bounds .upper  =  [tO  tf]  ; 
knots . numNodes  =  [150]; 


7. - 

7.  Set  variable  bounds 

% - 

bounds . lower . states ( : , 1) 
bounds .upper . states ( : , 1) 
bounds . lower . states ( : , 2) 
bounds . upper . states ( : , 2) 
bounds . lower . states ( : , 3) 
bounds . upper . states ( : , 3) 


[wO  hO  rO] ’ ; 

[wO  hO  rO] ’ ; 
[wlb  hlb  rib] ’ ; 
[wub  hub  rub] ’ ; 
[wlb  hlb  rib] ’ ; 
[wub  hub  rub] ’ ; 


bounds . lower . controls  =  -inf  *  [111]’; 
bounds .upper . controls  =  inf  *  [111]’; 


% - 

7.  Set  constraint  bounds 

7. - 

bounds . lower .path  =  [0] ’ ; 

bounds .upper .path  =  [le2]  ’ ; 

bounds . lower . events  =0*  [111111]’; 

bounds .upper . events  =0*  [111111]’; 


7. - 

7.  Provide  a  guess 

7. - 

load  iss_law_guess ; 


7. - 

7.  Call  DIDO 

7. - 

[cost,  primal, dual]  =  dido(iss,  knots,  bounds,  guess); 


The  above  main  file  is  similar  to  the  one  in  Section  B.l  because  the  control  in  this 
case  is  a  continuous  variable.  The  main  distinction  between  this  problem  set-up  and 
previous  ones  shows  up  in  the  differential  equation.  The  function  which  implements 
the  differential  equation  is  shown  below. 
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function  [R]  =  issDynamics (primal) 

7. 

°/„  ISS  momentum  control  dynamics 

l 

°/„  Jw’  =  Tgg(r)  -  w  x  (Jw+h)  -  u 

°/„  h’  =  kl*J*w_e  +  k2*J*r_e 

°/„  r’  =  0.5  *  (  rrJ  +  I  +  skew(r))(w  -  w(r)) 

°/„  global  variables 

global  ws  hs  rs  us  kl  k2  J  iJ  w_orb 

°/„  initialize  variables 
w  =  primal . states (1 : 3 ,: )  *  ws; 
h  =  primal . states (4 : 6 ,: )  *  hs; 
r  =  primal . states (7 : 9 ,: )  *  rs; 
u  =  primal . controls  *  us; 

N  =  length (w( 1 , : ) ) ; 

R  =  zeros(9,N) ; 

for  i  =  1:N 

%  compute  auxiliary  values 
C  =  eye (3)  +  2/(l+r ( : , i) > *r ( : , i) ) *  •  •  • 

(skew(r ( : , i) ) *skew(r ( : , i) )  -  skew(r( : ,i))) ; 

C2  =  C(: ,2) ; 

C3  =  C(: ,3) ; 

wp  =  -w_orb  *  C2; 

Tgg  =  3*w_orb~2*(cross(C3, (J*C3))) ; 

U  =  eye(3)  +  2/(l+u( : , i) J *u( : , i) )* . . . 

(skew(u( : , i) ) *skew(u( : , i) )  -  skew(u( : , i) ) ) ; 

Ct  =  C  *  U’; 

wt  =  w(:,i)  +  w_orb  *  U(:,2); 

pt  =0.5*  [Ct  (2,3)  -Ct  (3,2)  ;  Ct  (3 , 1)  -Ct  (1 ,3)  ;  Ct  (1 , 2)  -Ct  (2,1)]  /... 
sqrt (1+Ct (1 , 1) +Ct (2,2) +Ct (3,3)) ; 

%  compute  differential  constraint 

R(:,i)  =  primal . statedots (:, i)  -  [  iJ*(Tgg  -  ... 

cross (w( : ,i) , J*w( : ,i)+h(: ,i))  -  u(:,i))/ws; 

(kl*J*wt  +  k2*J*pt)/hs; 

(0.5*(r(: ,i)*r(: ,i) ’  +  eye(3)  +  skew(r ( : , i) ) ) * (w( : , i) -wp) ) /rs] ; 


end 


